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Abstract 

Disentangling the processes leading populations to extinction is a major topic in ecology and 
conservation biology. The difficulty to find a mate in many species is one of these processes. 
Here, we investigate the impact of self-incompatibility in flowering plants, where several inter- 
compatible classes of individuals exist but individuals of the same class cannot mate. We model 
pollen limitation through different relationships between mate availability and fertilization success. 
After deriving a general stochastic model, we focus on the simple case of distylous plant species 
where only two classes of individuals exist. We first study the dynamics of such a species in a 
large population limit and then, we look for an approximation of the extinction probability in small 
populations. This leads us to consider inhomogeneous random walks on the positive quadrant. 
We compare the dynamics of distylous species to self-fertile species with and without inbreeding 
depression, to obtain the conditions under which self-incompatible species could be less sensitive 
to extinction while they can suffer more pollen limitation. 

Keywords: birth and death process, ODE approximations, inhomogeneous random walk on the pos- 
itive quadrant, inbreeding depression, extinction probability, mating system, distyly. 
AMS Codes: 92D40, 92D25, 60J80. 



1 Introduction 



We are interested in modeling the specific sexual mating system of a plant population, and especially 
in highlighting several phenomena which can affect its dynamics: stochasticity in the demography, 
pollen limitation, boundary effects. Firs t, the fate of small populations depends on stochastic pro- 
cesses such as demographic stochasticity ( Lande . 19981 ). which refers to a variance due to randomness 
in the death and reprod uction events. Sec ond, when populations are small or at low density, there 
can happen Allee effects ( Allee et all . 19491 ). i.e. when a positive relationship between the size of the 
population and the per capita growth rate appears. In sexual species, the birth r ate may depend 
on th e availability in mating partner, what is called the "mate finding Allee effect" ( Gascoigne et all . 
20091 ). For example, in species with s eparate sexes, the mean reproductive rate of a population may b e 
low because the opposite sex is rare (jEngen et a 1, 2003; ISaether et all . 120041 : iBessa-Gomes et a l l2004h . 
These biases can be expected to be higher in small populations because of stochasticity. Finally, at 
the extreme, the mate finding Allee effect can result in boundary conditions where compatible mates 
disappear, leading inevitably the population to extinction. 

Many hermaphroditic species of Angiosperms (flowering plants) have mating systems which allow 
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fertilization only am ong specific classes of in dividuals: between long-styled and short-styled plants 
in distylous species ( Barrett and Shore . 20081 ). or amon g different matin g types in self- incompatible 
species (more than 50 percent of Ansgiosperm families, Igic et al ( 20081 )). Self- incompatibility (SI) 
avoids selfing and mating between close individuals: SI prevents reproduction between mates sharing 
identical S-locus alleles (the allele which determines whether individuals are compatible, by coding 
for some proteins carried as identifiers by the pollen and ovules), and especially self-fertilization. It 
was hypothesized that this mating system evolved to avoid inbreeding depression , i.e. the decrease 
in fitness when mating occurs between kin individuals ( Porcher and Landel . 2005al lbl). In this system, 
individuals carrying rare S-locus alleles have access to more possible mates than individuals carrying 
comm on S-alleles, which generates negative frequency-dependent selection on the S-locus (jWrightl . 
193flh . 



In Angiosperms, mate finding is moreover mainly passive since it depends on pollen vectors such as 
insects or wind. This can result in "pollen limitation", when a given plant does not receive enough 
pollen to fertilize all the ovules it produces. Pollen limitation has been found in many species and 
can have many evolutionary and ecological consequences, from the evolution o f mating system to the 
increase of extinction probabilities in small populations ksWnetaJ, H). 

If a class of individ- 
uals goes extinct, other classes have less mating opportunities leading to a lower mean reproduction 
rate, or even can not mate anymore, leading to extinction. Pollen limitation and matin g systems are 



i ntera cting phenomena which can have a strong impact on the fate of populations (see iLeducq et al 



( 201ol ) and references therein). Shortly, pollen limitation can be increased in strict allogamous species 
and the probability of disappearance of a class of individuals can be increased by pollen limitation. It 
has been shown indeed th at pollen limitation is higher in outcrossing species than in selfing species 
( Larson and Barrett . 2000l ). 

However, few theoretical investigations have been done so far to measure the impact of pollen lim- 
itation, mating systems and their interactions on the extinction and establishment of populations. 
It is yet a central question since many Angiosperm families are strictly allogamous because of Si - 
Thre e demographic models have been investig ated specifically for a given species (IKirchner et all . 
2006 : Wagenius et al . 2007 : Hoebee et al 20081 ). which thus lack generality. Levin et al ( 2009 ) con- 
sidered as us a gene ra l model for which t h ey performed indiv idual based simulations. Moreover, in 



Kirchner et a 1 (|2006h : iHoebee et all (|2008l ): iLevin et a 1 (120091 1. the supposed SI systems were either 



gametophytic (only the content of the gametes is expressed) or a caricatural sporophytic system with 
only codominance between S-alleles. In species with sporophytic SI system (SSI), the mating pheno- 
types of pollen and pistils are determined by the diplo id parenta l geno types at the S-locus and hence 
dominance interactions are p ossible among S-alleles ( Bateman . 19521 ). which can be very complex 
( Castric and Vekemansl . 2004] ) . The results from the previous studies can therefore not be generalized 
to species with SSI, and especially not to distylous species, which is a particular case of SSI with only 
two alleles and consequently only two classes of individuals. Finally, the different processes affect- 
ing the extinction probabilities of SI populations (pollen limitation, demographic stochasticity and 
boundary effect) can not be disentangled in these previous investigations, which is our purpose here. 



Our goals are, first to develop a general model to describe the dynamics of a plant population with 
SSI, with and without pollen limitation and second to use this model to investigate the relative impact 
of pollen limitation and demographic stochasticity on the fate of populations in the particular case of 
a distylous species. We begin with a general stochastic individual based model, in continuous time, for 
SSI (Section [2]). We explicitly model the genetic determinism of the SI phenotype and compute the 
reproduction rate of each possible genotype. We assume different relationships between the compatible 
mate availability and the reproductive success, which reflect different models of pollen limitation. 
Second, we focus on the simple case of distylous species. We analyze the dynamics of distylous species 
assuming a large population and using approximations by ordinary differential equations (ODEs). We 
exhibit different behaviors depending on the relationships between the birth and death parameters. We 
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also compare these behaviors with self-fertile species, with or without inbreeding depression and pollen 
limitation (Section [3]). Third, we consider extinction in small distylous populations (Section 2]). This 
leads us to study inhomogeneous random walks on the positive quadrant. We use coupling arguments 
to show that the behaviors of the random walks, namely whether they are subcritical, critical or 
supercritical, depend on the same criteria as for large populations. This also provides estimates for 
the extinction probabilities. Finally, individual based simulations are performed (Section [5|). 



2 Microscopic modelling 

We first describ e the individual dynami cs. Then, we precise the different models of reproduction 
rates. Following IChampagnat et a f tod ), we propose Stochastic Differential Equations (SDEs) that 
mathematically describe the random evolution in time of the population and its large population 
approximation using ODEs. 



2.1 Description of the dynamics at the individual level 

Self-incompatibility genotype We assume that SSI is controlled by a gene at a single locus, 
called the S-locus, where re alleles segregate. These alleles are numbered from 1 to re. We assume that 
individuals are diploid and hermaphroditic. Each individual genotype is of the form G = {S^S* 2 }, 
where S 1 and S 2 are two alleles in [1, re] = {1, . . . , re}. Since the order of the alleles is not important, 
for u, v G [l,re], {u, v} and {v,u} are the same genotype. 

We denote by E = {g = {u, v}, u,v £ [l,re]} the space of the different possible genotypes. The set E 
is finite, with Card(£) = n(n + l)/2. 



Individual based model The population at t is given by a point measure on E 

N t 

Z t {dg) = ^26 {sHils2m (dg) = Yl N t V '5{u>,v>}{dg) (2.1) 

i=l g'={u',v'}&E 

where Nt is the number of individuals alive at time t, Nt' v ' is the number of individuals with genotype 
{u',v'} G E. Each individual i is represented by a Dirac mass weighting its genotype {S 1 (i) , S 2 (i)} . 
Point measures that are considered have finite mass. We denote by Mf{E) the set of finite measures 
on E. Since E is a finite space, the weak and vague topologies, and the topology associated with the 
total variation norm on Mp(E) all coincide. 

Integrating the measure (|2.ip with respect to chosen real functions / on E provides summaries of 
the population. We denote (Z t J) = f E f(g)Z t (dg) = £^ f({S l (i), S 2 (i)}). 

If we choose for instance / = 1 then, (Zt, 1) = Nt- If we choose / = TLi u > v n for {u',v'} G E, then 
i z t, ![«',«'}) = Nt v ■ 

Self-incompatibility phenotypes SSI phenotype is the production of recognition proteins at the 
surface of pollen and stigmas (stigmas contain ovules and are the plant's structure on which pollen 
is deposited). These proteins can be of re different types (or specificities) depending on the plant's 
genotype and on the dominance relationships between the re alleles. For u G [l,re], let us denote by 
e u G M. n the vector with all components equal to zero except the component u which is equal to 1. For 
an individual of genotype {re, v}, the stigmas and pollen it produces are said to be of type {u, v}, and 
have respectively the phenotypes &s(eu + &v) and &p(e u + e v ) in {0, e u , e v , e u + e v }. The maps $p 
and $5 encode the expression of the genotype into the phenotype. For instance, if &p(e u + e v ) = e u , 
then u is strictly dominant over v in pollen (proteins produced at the surface of pollen are of the 
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single type u). If <&p(e u + e v ) = e u + e v , then u and v are codominant in pollen (proteins produced at 
the surface of pollen are of both types u and v). If &p(e u + e v ) = then the pollen can fertilize any 
possible ovule in the population. A cross is incompatible if pollen and stigmas share proteins of the 
same type. In other words, a stigma {u,v} and a pollen {u',v'} can cross if and only if 

$s{e u + e v ) ■ <S> P {e u , + e v >) = 0, (2.2) 

where x ■ y is the scalar product of x,y G W 1 . Notice that an individual can self-fertilize if 3>s(e u + 
e v ) • $p(e u + e v ) = 0. Let: 

8^ = {{u', v'} G E, $s(e u + e v ) • $ P (e u , +e*) = 0} (2.3) 

denote the set of genotypes compatible with stigmas {u, v}. The size of the population producing 
pollen compatible with stigmas {u, v} is hence: 

Wt V = N t V '=! U S {e u+ e v ).* p{ e u , + e v ,)=oZt{d{u' ,V'}). (2.4) 

{u',v'}eef Je 

Notice that the set 0™ of genotypes compatible with pollen {u, v} is not necessarily since we 
may have $s(e« + e v ) ■ ®p(e u > + e v >) ^ &p(e u + e v ) ■ ®s(e u ' +e v >). 



Mating probabilities Each plant receives pollen from the rest of the population. We assume that 
the quantity of pollen of type {u, v } G E received by a plant is equal to the proportion of plants 
of genotype {u, v} in the population. The probability that the ovule of an individual with genotype 
{u',v'} is fertilized by a pollen from an individual with genotype {n, v} is denoted by p uv (u,v). 
This probability takes into account the quantity of pollen received by the pistil and the compatibility 
conditions (|2.2|) . This is detailed in Section I2T21 



Segregation When the genotypes of the plants having produced the ovule and pollen are respectively 
{u, v} and {u',v'}, then the new individual is of genotype {n, u'}, {u,v'}, {v,u'} or {v,v'} with 
probability 1/4. 

Reproduction rates Plants {u, v} produce ovules in continuous time with an individual rate f > 0. 
Once produced, the ovules may be fertilized or not depending on the quantity and compatibility of 
the pollen arriving on the stigmas, thus providing seeds. We define as R(N t , Nt) the individual 
reproduction rate of an individual of genotype {u,v}, which is the individual rate of production of 
seeds. This rate is bounded by f and can depend on the size of the compatible population or be a 
constant: it is the product of r with the fertilization probability 

R{NT,N t ) = f Pt V (u',v>). 

{u',v'}&E 

The functional forms of R{.) are further discussed in Section \2.2\ 



Death Each plant dies with the constant rate d > 0. This death rate is kept simple and the 
complexity of the model is put on the reproduction system. 



2.2 Functional forms of mating probabilities and reproduction rates 

We describe the five typ es of reproduc tion considered in this paper. Three of them describ e SI 
systems: Wright 's model (jWrightl Il93fll ). the fecundity selection model (|Vekemans et all . [l998) an d 
the dependence model, which we introduce. In these models, ovules are produced at a constant rate, 
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and given reproduction, the partner is chosen uniformly among compatible individuals. The differences 
are based on the choice of functional forms for mating probabilities. The resulting reproduction rates 
for these three models are represented in Fig. [TJ Two other models are also introduced for comparison: 
the self-compatible cases with and without inbreeding depression. 

Boundary conditions, i.e. discontinuities of the reproduction rate that occur when a genotype in the 
population has no more mate, are observed in the Wright's and dependence models. For the fecundity 
selection model, even if there is SI, the reproduction rate is proportional to the size of the compatible 
population and vanishes continuously at the boundary. In self-compatible cases, there is no boundary 
effect because an individual can mate with any other individuals, even itself. 

In Wright's model, there is no pollen limitation. In the fecundity selection and the dependence models, 
pollen limitation is introduced through the mating probabilities. 



2.2.1 Model 1: Wright's model 

First, following Wright ( 19391 ). each individual with genotype {u, v} G E produces ovules of type 



{u, v}, one at a time, at the constant f. An ovule of type {u,v} is fertilized by a pollen produced by 
an individual of genotype {u', v'} with a probability that depends on the frequencies in the compatible 
population: 



,,/„,/ 



Pt v (u ,v ) = -=mr'&{u>,v'}eef>- ( 2 -5) 
^ t 

Notice that because of (|2.4p . pf v (u',v') = when N t = 0. The fertilization probability is 1 if there 
is at least one compatible individual in the population, otherwise. The individual reproduction rate 
is R(N t ,Nt) = ft-jfuv This rate does not depend on the quantity of compatible pollen a plant 
receives. 

2.2.2 Model 2: Dependence model 

Pollen limitation and the self-incompatibility may be modeled by considering mating probabilities of 
the form: 

r(N u '°) 

P r(u',v')=pr(u',v')^l (2.6) 
where pf v (u',v') is defined in (|2.5p and where, for positive constants a and /3: 

aN 

••W = fj^R- (2.7) 

In (|2.6p . the probability of finding a mate depends on the number of compatible individuals, and given 
reproduction, the partner is chosen uniformly at random among the latter. Thus, the fertilization 
probability is ij^^ >0 r(N t )/f and the individual reproduction rate is R(N t ,N t ) = r(N t )^jf uv >0 - 
This model will be termed "dependence model" in the sequel. It is a generalization of Wright's model 
with a fertilization probability that depends on N^ v when reproduction is allowed. This is pollen 
limitation: among the r ovules produced in a time unit, only R{N t ,Nt) end up in producing a new 
individual. For a large compatible population, the rate tends to its supremum value: 

e aN 

lim r— rrlLyx) = r. 

n^+oo j3 + e aN 

Wright's model can be viewed as the limiting case of (|2.7|) when a — > +oo, that is when there is no 
pollen limitation. 
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Wright's model 
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Qependence model / 

^ Fecundity selection model 



TT 



JVt-1 



Figure 1: Relationships between the individual reproductive rates of genotype {u, v} and the number 
of compatible individuals N ' for the three models of mating we assumed: Thin line: The Wright's 
model; Dashed line: The fecundity selection model; Thick dashed line: The dependence model. 



Notice also that this model exhibits boundary effects since 

e aN = 

lim f— w lljv>o = ~ - > 0. 

n^o (3 + e aN 1 + 

Because we assume that a single plant produces a lot of pollen, even in cases where we consider pollen 
limitation, there is a discontinuity at the set {N = 0} when r(0) > 0. As soon as there exists a 
small number of compatible plants, all individuals start producing offspring at a positive rate. 



2.2.3 Model 3: Fecundity selection model 



In the fecundity selection model ( Vekemans et al . 19981 ). the mating probability between an ovule 



produced by an individual of genotype {u, v} and a pollen produced by {u',v'} is 



pT{u\v')=p U t V {u',v') 



N 



N f 



^7-ll{«>'}eeF- 



(2- 



If the ovule chooses an incompatible pollen, then it is not fertilized and lost. Under the fecundity se- 
lection model, the fertilization probability of an ovule {u, v} is N t /Nt and the reproduction rate of an 
individual {u, v} is R(N t , iVj) = fN t /Nt. The individual reproduction rate is directly proportional 
to the frequency of compatible individuals, which thus promotes pollen limitation. 



2.2.4 Rates of reproduction and gamete production 

In the models of Sections 12.2.11 12.2.21 and 12.2.31 the rate at which pollen produced by individuals of 
genotype {u, v} fertilizes ovules in the population Z t is: 

£ rN^'p uW (u,v) (2.9) 
{u',v'}eE 

where we recall that p wv ' (u, v) is the mating probability of an ovule produced by an individual {«', v'} 
by a pollen produced by an individual {n, v}. Under Wright's model p u v (u,v) =p uv (u, v). 
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Hence, offspring with genotype {u, v} are produced with the rate: 



1 / \ - 



+2 1 2^ + E ^V™' 



u,u 



+rl u -- 



+ (Nr P uu (v,v)+Nrp vv (u, U ) 

E N? u 'p uu '(u,v') 

u',v'€ll,n}\{u} 

+\{ E ^rV u '(«,«) + E iv rp uu ( 

+N? u p uu {u,u) 



U, V 



(2.10) 



The formula (|2.10p does not simplify in general because of the possible asymmetry of dominance 
relationships between alleles in pollen and stigmas. We distinguish whether the allele u comes from the 
pollen or the pistil, and similarly for v. We separated the terms according to the parents' homozygosity 
or heterozygosity. 

The rate (|2.1U|) is related to the number of seeds of each genotype {u, v} that is produced at time 
t. There is a mass dependence with respect to the parent who brings the ovules and a frequency 
dependence with respect to the one who brings the pollen. The case of infinite populations with the 
Wright's and the fecundity selection model is considered in lBilliard et all (|2007l ). 

We end this section, with two alternative models without SI, for future comparisons with the three 
models introduced in Sections 12.2. Ill2~2. 31 Model 4 describes a self-compatible case without any penal- 
ization of self-fertilization, whereas Model 5 introduces inbreeding depression, th e phenomenon which 
descri bes the decrease in fitness when mating occurs between kin individuals (jPorcher and Lande . 
2005al Jbh. 



2.2.5 Model 4: Self-compatibility without pollen limitation nor inbreeding depression 

To carry out comparisons, it is natural to investigate the case where there is no SI, no pollen limitation 
and no inbreeding depression. Then, when the reproduction r is constant, the rate of production of 
offspring with genotype {u, v} becomes: 

r av (z t ) =|- (1 e N r' + nt) (I E N ™' + N t v ) iiu ^ v 

^ U{Zt) = w t (l ^ Nr ' + Nr T iiu = v - (211) 

In this rate, we recognize the number of couples that can be constituted with one parent of allele u and 
one parent of allele v, YZ'=i N t u E"'=i N ™ , divided by the total size of the population because of the 
frequency dependence, and with a correction for homozygous individuals. Notice that the genotypic 
frequency appears naturally since here, pf v (u',v') = iV" ' v /Nt is the proportion of individual {u',v'} 
in the population. 
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2.2.6 Model 5: Self-compatibility with pollen limitation and inbreeding depression 



A detailed study of inbreeding depression would be interesting and deserves a paper for itself. Our 
purpose here is to provide a schematic point of comparison for SSIs. For the sake of simplicity, we 
assume here that inbreeding depression is suffered by offspring produced by self-fertilization only. 
Among the ovules produced by the plant, a fraction s G [0,1] are self- fertilized, i.e. are fertilized 
by the pollen of the plant that has produced them. Among these, a fraction 5 G [0, 1] is not viable. 
The rate of production of new individuals through other matings is of the form (|2.7|) except that the 
maximal rate is now (1 — s)f which corresponds to the ovules that have not been self- fertilized. We 
hence consider the following rate in the case of pollen limitation: 



R(N t ) = (1 - S)sr + (1 - s)r 



e a(JV t -l) 

P + e «(ft-i) 



iV t >l- 



(2.12) 



In the sequel, this model will be termed by Self-Compatibility with Inbreeding Depression model 
(SCID). In the presence of self-fertilization, the effect of pollen limitation is lower: the reproduction 
rate of individuals is less dependent on the quantity of pollen received from the other individuals. 
With the rate (|2.12p . the rate of production of offspring with genotype {u,v}, at time t is: 



r« v {Z t ) =2(1 - a) 



MNt-l) 



N t 8 + e< N t- 1 ) V2 ^ 



if u 7^ v 



(2.13) 



^(Zt) =(1 - s) 



Ntp + ePt"*- 1 ) V2 

u'fu 



1 



N t U> + N t U ) + (1 - $)srN? u if 



U = V. 



2.3 SDE description and their ODE approximations 



Following iFournier and Meleard! (|2004h . we can express the evolution of the population {Zt)teM + by 
mean of a SDE driven by a Poisson point measure. This equation is given in Appendix IA.2I and 
corresponds to the mathematical formulation of the individual-based algorithm of Section 15.11 From 



this, we derive the evolution of N^ v : 



N uv =N uv 



+ f (r uv (Z s ) -dx N™)ds + M™ 
Jo 



(2.14) 



where M uv is a square integrable martingale starting at that can be written explicitly in term of 
the Poisson point measure and with bracket: 



{M l 



I \r uv (Z s ) + dxNr)ds. 
Jo 



(2.15) 



Heuristically, we can interpret the martingale as a noise term corresponding to the stochasticity and 
whose variance is given by (|2.15p . The integral term in (|2.14p gives the growth rate of the population, 
as for usual ODEs of popula tion dynamics, which are mo r e usual in the biolog i cal lit erature (see f|2.20|) 
in the sequel). We refer to Ikeda and Watanabel ( 19891 ): Joffe and Metivier ( 19861 ) for introductions 
to such SDEs. 

These SDEs are linked in large populations with classical ODEs. In similar cases with such nonlinear 
dynamics, we know that the se ODEs arise as deterministic approximations of the SDEs under a certain 
large population limit (see Champagnat et al . 20061 : Ethier and Kurtz , 19861 : Fournier and Meleard . 
2004 : Tran . 20081 ). Apart from providing the mathematical frame under which the SDEs can be 
approximated by ODEs, limit theorems may be interesting for statistical estimation, in particular for 



linking individual-base d statistics with parameter estimates for the ODE (e.g. Blum and Tran . 2010l : 
Clemencon et~all2008h . 



For the limit theorem, we consider a sequence of populations with initial sizes of order K (Point (/ — i) 
of Definition 12. ip . To counterbalance this effect, we renormalize the size of the individuals in 1/K 
(Point II) so that the population mass remains of the order of a constant. The renormalization of the 
reproduction rate (/ — ii) is a kind of localization factor: if the size of the population is large, only 
the neighborhood of a given individual will affect its reproduction rate. 

Definition 2.1. I) We consider the following sequence of processes (Z k )k&i* , where N* = N \ {0}. 

(i) Let (Z^)KeN* be initial conditions such that there exists a deterministic finite measure £o G 
A4f(E) such that: 

lim Z o(? g ) = £ (dg) and 3e > 0, sup e((^M 2+e> ) < +oo, (2.16) 

K-*+oo K KeN* " A ' 

where = (Z^ , 1) is the size of the population described by Z^ . 

(ii) The rate of production of a given genotype g, r 9 (.) is replaced by r 9,K (.) = r 9 (./K). The death 
rate d is unchanged. 

II) We also introduce the sequence of renormalized processes (Z^)k^n* ■' 

G N*, yt G R+, Zf\dg) = ^Z t K (dg). (2.17) 

We define by = /K the mass of z[ K \ that is the renormalized size of the population. ■ 

We use exponents (K) with brackets for the renormalized quantities and exponents K without 
brackets for the non-renormalized ones . The moment condition in (|2.16p is technical and ensures that 
(for a proof, see Fournier and Meleard . 20041 ) : 

VT G R+, sup E( sup {N^ K) f) < +oo. (2.18) 

KeN* te[o,T] 

The ODEs are obtained when K — > +oo. The corresponding limit theorem is stated in the next 
proposition and proved in appendix. 

Proposition 2.2. The sequence of renormalized processes (Z^)KeN* converges uniformly, when 
K — > +oo to the process in C(M+, M.p{E)) such that: 

Z t {dg)= Y, < 6 {u,v}( d 9)> (2-19) 

g={u,v}£E 

where for every {u, v} G E, nf v is the continuous number of plants of genotype {u, v} at time t and 
satisfies the following ODE: 

dn uv 

-^ r =r uv {it)-dnT (2.20) 
where r uv (^t) is obtained from H2.10\) by replacing all the N uv 's by n uv 's. 
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2.4 Particular case of distylous flowers 



The most simple case of SSI is distyly where n = 2 and $p = $5 = We detail this case, which will 
be studied carefully in all the rest of the paper. We will see that the two alleles in this case can not be 
codominant, which corresponds to what happens in natural populations. This leads us to introduce a 
random walk on the positive quadrant, which is central in the sequel. 

The possible genotypes are {1, 1}, {1,2}, {2,2}. The population is described with N , iV 12 and 
iV 22 . Two alleles 1 and 2 are codominant if: 

$(2ei) = ei, $(2e 2 ) = e 2 , $(ei + e 2 ) = ei + e 2 . (2.21) 

The allele 2 is dominant and 1 is recessive if: 

$(2ei) = ei, $(2e 2 ) = e 2 , $( ei + e 2 ) = e 2 . (2.22) 

Proposition 2.3. In the case of codominance \2.21\) , there is almost sure extinction of the population. 

Proof. Since <&(ei + e 2 ) • $(ei) = $(ei + e 2 ) • <5(e 2 ) = 1 and $(ei + e 2 ) • $(ei + e 2 ) = 2, heterozygous 
individuals {1,2} are unable to mate with any individual. The only possible match is between indi- 
viduals {1, 1} and {2, 2}, but this produces individuals {1, 2} which have no compatible mate. Hence, 
there is at most two generations of individuals. Since the death rate is a positive constant d > 0, 
extinction happens in finite time almost surely. ■ 

In the sequel, we will therefore assume that alleles 1 and 2 are not codominant. Let us consider 
the case where 2 is dominant and 1 recessive (|2.22|) . the symmetric case being treated in the same 
manner. In this case, {1,1} can mate with {1,2} and {2,2} whereas the latter can only mate with 
{1,1}. The sizes of the respective compatible populations (|2.4|) are: 

Af=^ 12 + A? 2 , N? = Nl\ Nf = N}\ (2.23) 

First, since none of these matings produces offsprings of genotype {2, 2}, this genotype is only present 
in the first generation and hence disappears in finite time almost surely. For the sake of simplicity, we 
assume that the initial condition is only made of individuals of genotypes {1, 1} and {1,2}. In this 
case, the genotype {2,2} never appears and A^ 12 = N t — Nj Individuals of genotype {1, 1} can only 
reproduce with {1,2} and reciprocally. As soon as one of the genotypes {1,1} or {1,2} disappears, 
the whole population is doomed to extinction since the remaining genotype can not reproduce any 
more. We are led to study random walks on the positive quadrant, absorbed on the axis {x = 0} and 

{y = o}. 

In the absence of {2, 2}, p (1, 2) = p (1, 1) = 1. The random walk (|2.14p becomes: 

+ J Q (R(N?, N S )N? + R(N'\N S )N?) -dx ds 
N} 2 =N^ 2 + M} 2 (2.24) 
+ Io (li^s^^N^+RiN^ 1 ,^)^ 2 ) -dx N} 2 ^j ds. 

A particular significance is given to the escape time from the first quadrant: 

r = inf{t G R+, A^ 11 = or A^ 12 = 0}. (2.25) 

On the set {r < +00}, the population goes extinct in finite time almost surely, whereas on the set 
{t = +00} the population survives forever. 
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3 Large populations of distylous flowers 



We consider the Models 1-3 of Section HT21 for the distylous populations of Section [23] and address the 
questions of determining when SSI are advantageous with respect to self-fertilization and inbreeding 
depression. The study is carried in the case where the population is large and the ODEs studied in 
Section [2T3l are used. The case of small population is tackled in Section HI 

In Sections 13.11 to I3.3( we study the stationary solutions and their stabilities in the case of Models 
1 to 3 (Sections 12.2.11 to I2.2.3p . Threshold phenomena and asymmetries showing Allee effects are 
highlighted. In Section \'dA\ comparisons between a SSI and a system with self- fertility and inbreeding 
depression are studied. 

For Wright's (Model 1) and dependence (Model 2) models, the ODEs are: 



n? J ~ W 2 ) Jo V KrK^^ + r^K^V^V^-dn 12 » ^ 

where r(.) = f for Wright's model and r(.) is defined in (|2.7|) for the dependence model. This gives 
with classical arguments on the regularity of solutions: 

^ =\{r{n] 2 )n? + r(n* VP) l n u >0 l n x 2>0 - dn? 
dn} 2 1 



H 



-(r-K 12 K n + r(n t n )nP)i 11>0 i 12>0 - dn] 2 , (3.1) 



dt 2 

with the initial conditions (ng , rip 2 ). The roles of n 11 and n 12 are symmetric. It is easy to see that 
when n 11 vanishes, the solutions remain on the boundary {n 11 = 0}. There is existence and uniqueness 
of the solutions on (R* + ) 2 = (R+ \ {0}) 2 , {0} x R + and R + x {0}, from which we deduce the existence 
of a unique solution for every initial condition (rig 1 , rig 2 ). 

3.1 Wright's model (Model 1) in a large population 

The system (|3.ip without the indicators in the right hand sides (r.h.s.) and with r(.) = f becomes: 
l i^=A( n t\ where A = ( % ~ d f i ) = L J - dl, (3.2) 



dt V nf / \ n t J V i h~ d ) 2 

where J is the square 2 x 2-matrix filled with ones and I is the identity 2x2 matrix. 

Proposition 3.1. There exists a unique solution for h3.1\) . which coincides with the solution of {2 
for every t G R + : 

n\ l J- (nf (e( f - d )* + e~ dt ) + n 12 {e^ - e~ dt 
n\ 2 =I(ni 1 (e( f - rf ) t - e~ dt ) + n 12 (e^ + e"*)). (3.3) 

and: 

lim e-^n} 1 = lim e^-^n] 2 = (nf + nJ 2 )/2. 

t— >+oo t— Y+OO 

Corollary 3.2. We recover a natural disjunction between f > d, f < d and f = d. The first case 
is the supercritical case: the population survives and grows to infinite size. The second case is the 
subcritical case: there is asymptotic extinction of the population. In the third case, the critical case, 
the population size remains constant and the solution converges to a non-trivial equilibrium. 
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Figure 2: Solutions {n\\n\ )t>o of \3.1\) in the (a): critical (r = 2, d = 2), (b): supercritical (r = 4, 
d = 3 J and (c): subcritical (f = 2, d = 3) cases. 



Remark 3.3. It is remarkable that in the large population approximation, for constant reproduction 
rates, the behavior is the same as for the case of compatible reproduction, in which the ODE for the 
population size is: dn/dt = (f — d)n. ■ 

Proof of Prop. \3.1[ Existence and uniqueness of the solutions of ([32]) in C°°(IR + ,M 2 ) hold by Cauchy- 
Lipschitz theorem. In every cases, the solutions of (|3.ip coincide with the solution of (|3.2p until one 
of the components equals zero. Let us denote by t° the time at which this happens and which may be 
infinite in case the trajectory of the solution of (|3.2[) does not intersect the horizontal or vertical axes. 
This implies existence and uniqueness until the time t° of the solution of (|3.ip . 

Let us solve (|3.2p . The matrix A admits r — d and — d as eigenvalues respectively associated with 
the eigenvectors (1/^/2,1/^/2) and (1/^2, -1/^2). For the system (pT2"j) . there is a unique solution 
for every initial condition (rig , nj 2 ) G (0, +oo) 2 such that t t-t (n\ , n\ ) is of class C°°. 



nf V J_ ( 1 1 \( Coe^-W 
nf ) ^2 V 1 - 1 J V C x e- d - 1 



(3.4) 



The constants Co and C\ are obtained from the initial condition. Solving (|3.4p in (Co, C\) for t = 0: 

C = ^(n 1 +n 2 ), C^^K-no 2 ). (3.5) 

Hence this provides (|3.3p . 

The question is whether (R?j_) 2 is positive invariant, i.e. whether the trajectories of (|3,3p remain 
in the positive quadrant. 



nl 2 - nl l 



n t n = o nJ 1 (e« + l)+nS 2 (e w -l)=0 O e w = ^ ° 2 . (3.6) 

n Q + n 

This equation has no positive solution in t since the right hand side is smaller than 1 (and possibly 
non positive). Similar computation holds if one solves n\ 2 = 0. Hence, the solutions of (|3.ip and (|3.2p 
coincide on 1R+ and to = +oo. 

The long time behavior is obtained by noticing that whatever the case, the dominant factor in 
is exp((f — d)t). ■ 
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3.2 Dependence model (Model 2) in a large population 



We now turn to the case of a variable individual reproduction rate r(.). First, we consider the general 
case of a bounded continuous nonnegative function r(.), and then we focus on the functional form 
given in (|2.7p . The ODEs limits are given in (|3,ip . Conclusions are summed up at the end of the 
subsection. 



Let us start with a bounded continuous nonnegative function r(.). The trivial solution (0,0) is a 
stationary solution. We look for non trivial stationary solutions. Because of the symmetry in n 11 and 
n 12 , if a nontrivial stationary solution exists, it satisfies n 11 '* = n 12 '* := n*. The latter value solves: 

r(n*)n* = dn* r(n*) = d. (3.7) 

The number of non trivial fixed points depends on the number of roots of 



Remark 3.4. Notice that the total population size at equilibrium is then 2n* , which is twice the size 
at equilibrium in absence of SI. Indeed, in the latter case, the size of the population at equilibrium 
satisfies flg. 7|) and is thus n* . ■ 



We now examine the stability of the stationary solutions. We r efer to IVerhulstl (|2000h for definit 
and developments on the theory of dynamical systems. 

Proposition 3.5. The trivial equilibrium (0,0) is: 

(i) a positive attractor if r(0) < d, 

(ii) a saddle point if r(0) > d. 



ions 



Proof. Because of the indicators in (|3.ip . we know that once one component has reached zero, it 
can not escape. We consider the stability of the trivial solution (0, 0) for the ODE (|3.ip without 
the indicators, as we know that before one of the components reac hes zero, these systems have the 
same solutions. We use the classical linearization methods {e.g. Verhulst . 2000l . Chap. 3). The 



linearization of the ODE without indicators around an equilibrium (re n ,n 12 ) leads us to consider the 
Jacobian matrix of the system at this point: 



. / 176 )IL J / \ll III. 

~ll 11 12\ _ f 2 ' 2 a 2 ' 2 — I (1 o\ 

J\n , n J- | r(n 12 ) | r'(n n )n 12 r'(n 12 )n n , r(n u ) , I ^' 8 >> 



r(n 12 ) j_ r'(n 11 )n 1 ' 2 ^ r'(n 12 )ra 1:L r(n 1:L ) 

r'(n n )n 12 r'jn 12 )^ 11 _, 

\ 2 + 2 2+2 

For the equilibrium (0,0), with the notation J and / introduced after (j3.2|) : 





:U<U»)- ( ^ (0 / r( T J )= r -^J-dI. (3.9) 



This matrix is the same as the matrix A in (|3.2p and its eigenvalues are r(0) — d and — d respectively 
associated with the eigenvectors (l/\/2, 1/^/2) and (1/V2, — 1/\/2). If r(0) < d, then both eigenvalues 
are negative and (0, 0) is a positive attractor for the system without indicators, and hence also for 
(|3.ip . If r(0) > d then there is a positive and a negative eigenvalue. In this case, (0, 0) is a saddle point 
for the system without indicators. This entails as in the proof of Prop I3TT1 that n\ l +nj 2 converges to 
+oo while rij 1 — n\ 2 converges to zero. Thus, in the neighborhood of zero, when starting from points 
of the positive quadrant, the solutions are the same as for the system without indicator: there is no 
extinction. ■ 

Proposition 3.6. For an equilibrium (n*,n*) with n* > 0: 

(i) Ifr'(n*) < 0, then the equilibrium is a positive attractor, 

(ii) if r'(n*) > 0, then the equilibrium is a saddle point. 
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Proof. For the point (n*,n*) the Jacobian matrix of f)3.8[) becomes: 



Z(n*,n*)=( r -^ + r 'K )n *)J-dL (3.10) 



This matrix is again of the same form as A introduced in (|3.2[) . Its eigenvalues are r{n*)-\-r'{n*)n* — d 
and — d, respectively associated with the eigenvectors (l/y/2, l/\/2) and (l/y/2, —l/y/2). As r(n*) = d, 
we can simplify the expression of the first eigenvalue: r(n*) + r'{n*)n* — d = r'(n*)n* , which is of the 
sign of r'{n*). As usual, the sign of the eigenvalues determines the nature of the equilibrium. ■ 



Let us now focus on the particular form (|2.7p for the reproduction rate. We compute the stationary 
solutions (n*,n*) of (|3.ip by starting from (j3.7j) : 



e an 1 / Bd \ 

'•(«*) = d * r 7^Tp = d & n ' = a l <V^)- (3 ' n) 

Of course, we see that the log is well defined if and only if f3 > and r > d. Moreover, (n* ,n*) belongs 
to the positive quadrant if and only if 

Ul >1 r<d{l + P). (3.12) 



r — d 

Since r(0) = r/(l +/?), we notice that the stability condition for the trivial equilibrium, r(0) < d, 
is equivalent to the condition (|3.12|) for the existence of a non-trivial stationary solution. Hence, if 
r(0) < d, (0,0) is a positive attractor and there is no other equilibrium, and if r(0) > d, (0,0) is a 
repulsive attractor. Indeed, (0, 0) is a saddle point for the system (|3.ip without the indicators. Since 
the stable manifold is the vectorial line of direction (1,-1) which intersects the positive quadrant only 
at (0,0), then for the system f)3. 1 j) . the equilibrium (0,0) is a negative attractor. 

Moreover, the equilibrium (n*,n*) is always a saddle point as: 



>(f3 + e an )-ae 2an raf3e c 
rW " r (/3 + e<™) 2 -(/3 + e <™ )2 > (•'■'•>» 



is always positive. The stable (resp. unstable) manifold is locally the affine line of direction (1, —1) 
(resp. the affine line of direction (1, 1)). 
In conclusion: 

• If f < d, then every trajectory converges to (0,0). 

• If d < f < d(l+/3), then there exists a non-trivial equilibrium that is a saddle point. Trajectories 
converge to (0,0) or lim t _;. +00 n\ 1 = lim^+oo n| 2 = +oo. 

• If f > d(l + /3), then there is no non-trivial equilibrium in the positive quadrant as the growth 
rate is too strong. (0,0) is a negative attractor and lim t _;. +00 n\ 1 = lim^+oo n\ 2 = +oo. 



Remark 3.7. Threshold phenomena appear: when d < f < d(l + /3), we see that contrarily to the case 
of absence of pollen limitation (Wright's model), there may be either survival or extinction. If the size 
of the population is too small (if the initial condition belongs to the attracting domain of (0,0)) then 
there is extinction. There is a threshold implying that survival is possible only for sufficiently large 
population. 

On Fig. 0, we see that the minimum size to ensure survival also depends on the composition of the 
population. It is smaller for population with 1/2 individuals of genotype {1,1} and 1/2 individuals of 
genotype {1,2}. If we draw the affine line of direction (1,-1) going through {n*,n*) on Fig. [3] (c), 
it can be seen for a given initial condition uq = raj 1 + raj 2 close to 2ra* that the trajectories either 
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Figure 3: Solutions (n^,nj 2 )t>o of (E2P / or a reproduction rate of the form {2.1^ with a = f3 = 1. 
fa): r = 2, d = 2>, (b): r = 7, d = 3. We have n* = —0.3 and this equilibrium is not in the 
positive quadrant. r(0) = 3.5 > d and (0,0) is a negative attractor. Every solution tends to infinity, 
(c): r = 3.1, d = 3. We obtain n* = 3.4. Since r(0) = 1.55 is smaller than d = 3, (0,0) is a 
positive attractor in this case. The solutions hence either converge to zero or to infinity, given the 
initial condition. It is also seen that at the neighborhood of 2n* , for similar initial population sizes, 
populations with higher symmetry are more likely to survive than very asymmetric populations. 



lead to +oo or when t — > +oo. Thus, the symmetry in mating partners matters for survival or 
extinction (Fig. (c)). To ensure survival whatever the initial condition, one needs a reproduction 
rate sufficiently large (f = d(l + /3)) to totally compensate the pollen limitation. 

Both phenomena can be interpreted as Allee effects: due to pollen limitation, the population growth 
rate is an increasing function of the compatible population size and thresholds appear when the function 
is negative at the beginning. ■ 

3.3 Fecundity selection model (Model 3) in a large population 

We finally study Model 3. The rate of production of genotypes {1, 1} is r , which does not present 



N t 

is now: 



any discontinuity on the boundaries {N^ 1 = 0} and {A^ 12 = 0} any more. The ODE approximation 



dn} 1 _ n\ x n} 2 , n dn\ 2 _ n} l n\ 2 12 

—t~ = r —\T — -v? ~ dn t\ — n~ = r -fT — ~v? ~ dn i > 3 - 14 ) 

dt n] 1 + n] 2 1 dt n] 1 + n\ 2 1 

which admits a unique solution in (M?j_) 2 . 

Proposition 3.8. As for Wright's model, we have a disjunction into three cases: 

(i) Subcritical case f < 2d: lim^+oo n\ l = lim t ^ +00 nj 2 = 0. 

(ii) Critical case f = 2d: the size n^ 1 + n\ 2 remains constant and {n\ l ,n\ 2 )t^M + converges to an 
equilibrium (n*,n*) £ (M^) 2 . 

(Hi) Supercritical case f > 2d: lim^+oon^ 1 = linn_ s>+00 n\ 2 = +oo. 

The additional coefficient 2 in the criteria comes from the fact that at equilibrium, the probability of 
mating success of an ovule is 1/2 contrary to the case of \3. 7\) where it is 1. 



Proof. By symmetry, a non-trivial equilibrim exists if and only if: 

rn* = 2dn* o f = 2d. (3.15) 
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To obtain the stability of the non-trivial equilibria, we compute the Jacobian matrix of the system 
(HS1D at (n n ,n 12 ): 

a(n n ,n 12 )=l ^7^,2 , V ™xV I (3-16) 

\ r \n™+n^ ) r ( n n n +n i2 ) ~d J 

For an equilibrium (n*,n*), = \ J — dl, and we are led to a discussion similar to (|3,10p . A 

difficulty arises at the trivial equilibrium (0, 0) since the limit of (n 12 /(n u + n 12 )) 2 at (0,0) in 5(0,0) 
is not defined. To remedy for this, we consider the ODEs satisfied by the total size nt = n\ 1 + n\ 2 
and by Zt = nj 1 x n\ 2 : 

dixit n} l n} 2 , dzt dn} 1 19 n n f 12 / \ , 

^^rip-*- l^ = ^r"' + "'^F = ( f - 2<i K (3 - 17) 

As a consequence, when f < 2d, there is asymptotic extinction of at least one of the genotypes {1, 1} 
or {1,2}, let us assume that it is n . Since: 

^ < fn\ l -dnf<e- dn] 2 (3.18) 

for all e > and for sufficiently large times, n 12 also converges to zero. 
If f > 2d, 

lim n} 1 x n} 2 = +oo. 

i->+oo 

One at least of the two subpopulation sizes {1, 1} or {1, 2} tends to infinity. Let us assume that it is 
n 11 and that n\ l > n\ 2 for sufficiently large t (implying that n\ jfit > 1/2). Writing: 

^>\n?-dnf, (3.19) 
we see that n 12 also necessarily tends to infinity. 

Finally, when f = 2d. If n] 1 < nj 2 , then nj 1 /2 < n] v n\ 2 jn t < nj 2 /2, and 

dn]l > f? - d]nV = and ^ < f? - ^n| 2 = 0. 



dt V2 / 1 (it V2 

Thus n 11 increases and n 12 decreases until they reach a state where n 11 = n 12 . ■ 

3.4 Comparison with self-compatibility with inbreeding depression (SCID) 

Now that our three models of interest have been considered, we wish to carry out comparisons with 
the SCID cases. Under which conditions do our models predict that SSI are advantageous (in some 
sense defined in the sequel) with respect to self-fertilization ? 

3.4.1 Self-compatible distylous population with inbreeding depression 

We consider here reproduction rates of the form (j2. 12[) . In the case of self-compatible populations 
with inbreeding depression, the size n t of the population in a large population limit follows the ODE: 



dn t 
dt 



(1 _ 5)sf + (1 - s)r JT -^ - d\ nt. (3.20) 



The trivial solution defines an equilibrium and: 
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Proposition 3.9. A non-trivial equilibrium to (|3,20p exists if and only if 

p[d- rs(l - 5)} > [r(l - 5s) - d] > 0. (3.21) 

(i) If IIS. 21]) is satisfied, the non-trivial equilibrium n* is given by: 

. 1 f P{d-(l-6)sf) ^ 

n = - log — r — — 3.22 

a \ r(l — os) — d ) 

This equilibrium is repulsive and is a locally positive attractor. If no < n* then the solution converges 
to zero, else, it grows to infinity. 

(ii) Else, there is no non-trivial equilibrium in the positive quadrant and there are three possibilities: 

• If r{\ — 5s) — d < 44> f(l — 5s) < d, then there is extinction (d is too large). 

• If fi[d—rs{\ - 5)] < O fs(l -5) > d, or if [r(l - 5s) - d] > f3[d — rs(l - 5)] > 0, then there 
is growth to infinity. 

Proof. A non-trivial equilibrium n* of (|3.20p satisfies necessarily: 

e an * _ , 1 / (3(d-(l-5)sr) 

[\ — s)r— = d — (1 — 5)sr 44> n = — log ' 



f3 + e an * a V(l - s)f - d+ (1 - 5)sfj 

which provides (|3.22|) if the log is well defined: 

(1 - 5)sr <d<(l- 5s)r. (3.23) 

By comparison arguments, we can prove that for d > (1 — 5s)r, every solution converges to as the 
square bracket in the r.h.s. of (|3.20p is strictly negative. For d < (1 — 5)sr, this bracket remains 
strictly positive and the solutions converge to +oo. Under the condition (|3.23p . n* is strictly positive 
if and only if 

(3(d - (1 - 5)sr) > (1 - 5s)r -d ^ (1 + (3)d > r (l - 5s + (3 s(l - 5)) 

P[d- rs(l - 5)] > [r(l - 5s) - d] . (3.24) 

When ()3.23p is satisfied, the brackets in (|3,24p are positive. Equation ()3.24p says whether the pa- 
rameter d in (|3.23p is closer to the lower bound of (|3.23[) {i.e. < f3(d — fs(l — 5)) < f(l — 5s) — d, 
implying growth of the population size to infinity) or of its upper bound (i.e. (|3.2ip . under which 
there exists a non-trivial equilibrium). Hence, we obtain (|3.2ip as sufficient and necessary condition 
for the existence of a non-trivial equilibrium in the positive quarter plane. 

With (|3.2ip . the situation is similar to (|3.1ip with d — (1 — 5)sr instead of d and (1 — s)r instead 
of r. Self-fertilization amounts to a reduction of the natural mortality since individuals can at least 
mate with themselves. However, this introduces a limitation of the maximal number of individuals 
produced by outcrossings. 

To study the stability of the equilibria and n*, we linearize the system and repeat the computation 
of (|3.9p and (|3.10p . The stability of n* depends on the sign of r(n*) + r'(n*)n* — d = r'(n*)n* which is 
here always positive by (|3,13p and by the remark of the preceeding paragraph. Hence, n* is a negative 
attractor. The trivial equilibrium is a positive attractor if and only if 

r(0) < d & (1 - 5)sr + < d & (l - 5s + /3s(l - 5))r < d(l + P). 

We recognize here the condition (|3.24p of existence of a non-trivial equilibrium, which leads to the 
announced result (i). In the case where (|3,2ip is not fulfilled, the behavior of the solutions is obtained 
by comparisons with simple ODEs. ■ 
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3.4.2 Comparison with the self-incompatible case 



It is natural to compare the self-compatible model of Section 13.4.11 with the SI distylous model of 
Section I3T2"! with the same reproduction rate (|2.7p . To carry the comparison, we introduce the following 
criteria: 

• We compare, when they exist, the population sizes at the non-trivial equilibria. Heuristically, 
these sizes provide the limit between the two extreme behaviors that are growth to infinity and 
extinction. When this equilibrium size is high, large populations are needed to avoid extinc- 
tion and the population is more "fragile". The corresponding model will be said to be less 
advantageous. 

• We can also compare the range of parameters r and d (for fixed s and 5) for which the population 
goes extinct. We will say that with respect to this second criterion, the more advantageaous 
conditions correspond to the smaller ranges of such parameters. 

We sum up the results of Sections 13.21 and 13.4.11 in Table [2j 



Case 


Case description 


Behavior of the population size 


Wright's model (Model 1) 


(a.l) 
(b.l) 
(c.l) 


f < d 
f = d 
f > d 


Convergence to 

Convergence to a non-trivial fixed point 
Divergence to +oo 


Dependence moc 


lei (Model 2) 


(a.2) 
(b.2) 
(c.2) 


r < d 

d<f<d(l + P) 
f > d(l + /3) 


Convergence to 
Existence of a saddle point 
Divergence to +oo 


Fecundity selection model (Model 3) 


(a.3) 
(b.3) 
(c.3) 


f <2d 
f = 2d 
f>2d 


Convergence to 

Convergence to a non-trivial fixed point 
Divergence to +oo 


Self-compatible model without in 


areeding depression (Model 4) 


(a.4) 
(b.4) 
(c.4) 


f < d 
f = d 
f > d 


Convergence to 

Convergence to a non-trivial fixed point 
Divergence to +oo 


Self-compatible with inbreeding depression (SCID) model (Model 5) 


(a.5) 
(b.5) 
(c.5) 


(1 - 5s)r < d 

p(d - fs(l - 5)) >f(l-5s)-d>0 
f > cZmin 1+/3s _ s5(1+/3 ) ) 


Convergence to 

Existence of a repulsive equilibrium 
Divergence to +oo 



Table 1: Summary of the behavior of the population size in the Sections \3.2\ and \3. 4 ■ 1\ depending on the respective 
values of f, d, s and 5. The condition for Case (c.5) is obtained by saying that we have divergence to +00 in 
the SCID model if (1 - S)sf > d or if (1 - S)sf < d < fs(l - 6) and f(l - Ss) - d > f3(d - fs(l - 6)) > 0. 

Models 1, 3 and 4 are similar. We have convergence to zero or divergence to infinity except in the 
particular cases when f = d or f = 2d. This shows that under fecundity selection, distylous species 
are more fragile than under Wright's model or self-compatibility without inbreeding depression. This 
is expected since there is no pollen limitation or self-incompatibility in the two latter models. Notice 
also that in absence of pollen limitation, the equilibria and critical, sub- and supercritical regions 
correspond in Models 1 and 4. 
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In the sequel, we compare Models 2 and 5, for which there exists a range of parameters for which 
extinction or divergence to infinity coincide, depending on the initial condition. In both cases, the 
pollen limitation is modeled similarly (see (|2.7p ). The difference relies on the penalization by SI in 
the first case and by inbreeding depression in the second case. This comparison provides conditions 
on the parameters under which distyly is advantageous on self-fertilization. 



Comparison with respect to the range of parameters The cases (a. 2), . . . , (c.2), (a. 5), . . . (c.5) 
are defined in Table [2j Let us consider the set of parameters for which the population goes extinct 
whatever its initial condition (Cases (a. 2) and (a. 5)). Since 5 > and s > 0, we see that this region is 
larger for the SCID model which in this respect appears as less avantadgeous than the SI model. The 
parameter (1 — 5s) in Case (a. 5) can indeed be interpreted as the proportion of seeds which survived, 
i.e. when excluding the fraction of non- viable seeds produced by self-fertilization. This term thus 
appears as an extra death parameter that is not present in the SI model. In our large population set- 
ting, this penalty is more important than the loss of partners that may face individuals in the SI model. 

We now turn to the Cases (c.2) and (c.5) where the population size diverges to infinity. The SI 
model is advantageous on the SCID model if and only if: 

(1 + ® ~ min ((1^7)? l + /3 S - + i(l + /3)) 
<=> j— >min(s(l-<5),l-<5) =s(l-S). 

This condition shows that for low self-fertilization efficiency, i.e. for large inbreeding depression, it 
is less advantageous than SI. The latter efficiency is expressed by comparing the fraction of viable 
seeds produced by self-fertilization to the initial fraction 1/(1 + = r(0)/f of seeds produced without 
self-fertilization when only an infinitesimal quantity of compatible individuals is present. 



Comparison with respect to the sizes of the population at equilibrium Let us now consider 
the cases when both (b.2) and (b.5) are satisfied. When there exists a non-trivial equilibrium, we have 
seen that the behavior of the solution is determined by its initial condition. The size of the equilibrium 
provides an idea of how many individuals are necessary to allow survival, even if in cases as in (|3.1ip . 
the symmetry of the initial condition may matter. Let us thus compare the sizes at equilibrium in 
(|3.1ip and (|3.22p . The size of the population at equilibrium in the SCID model is equal to the size of 
the population in the SI model when: 

a V (1 — os)r — a J a \r — dJ (1 — os)r — d \r — dJ 

Let us study: 

(3d-Pf{l-5)s 

fM = (r-d)-r5s ■ (3 ' 26) 

Notice that s = corresponds to the self-compatible case without inbreeding depression: /(0, 5) = 
(3d/ (r — d). The case s = 1 provides a model of compatible population with reproduction rate f (1 — 5) 
and we have f(l,5) = —f3. For any given 5 £ [0,1], f(-,5) is a rational function of s, defined on 
[0, 1] \ {so} with so := (r — d)/r5. The latter value belongs to [0, 1] if and only if 

d<r<— ^— . (3.27) 
l — o 



On the domain [0, 1] \ {so}: 

(3r(r(l -6)-d ) 

((r — d) — r5s) 
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which has the same sign as d — r(l — 5). 
In a nutshell, 

• If d < r(l — 5), the s i-> /(s, 5) is a decreasing continuous function on [0, 1], bounded above by 
/(0, 5) = /3d/ (r — d) < (/3d/ (f — d)) 2 (which is larger than 1 in Case (b.2)). In this case, there is 
no solution to (|3.25p . The size at equilibrium of the model with self-fertilization is always lower 
than its counterpart with self-incompatible reproduction: the cost of self-incompatibility may 
not be counterbalanced in this case. 

• If d > r(l — 5), then f(-,5) is an increasing function from [0,so) i n t° [/(0, £), +oo) and from 
(so, 1] to (— oo, — f}\. Thus, so is the highest fraction of self-fertilization that allows the existence 
of an equilibrium f|3.22|) and we will thus only consider s < sq. There exists a unique solution 
r G [0, so) to (|3.25p given by: 

(^l) 2 (r-d)-/3d 

r = 9 • 

For any s 6 (ro,so), the size of the equilibrium in the self-incompatible case is smaller than its 
counterpart with possible self-inbreeding. This means that if the fraction of offspring produced 
by self-fertilization is too high, then it is advantageous to switch to self-incompatibility. 



4 Small distylous populations 

We now consider cases where the deterministic approximation may not be taken and stick with random 
walks. In this section, we focus on the cases with constant rate of ovule production (Wright's model 
(Model 1), the fecundity selection model (Model 3) and their counterpart with possible self-fertilization 
(Model 4)) for which computations are tractable. Simulations for the other models are carried in 
Section [5j We are interested in the study of extinction probabilities. A difficulty comes from the fact 
that the random walks are inhomogeneous. After providing an equation for the extinction probabilities, 
we proceed by couplings to obtain approximations of these quantities. We construct processes defined 
on the same probability space as the original process (TV* , N{ } 2 )teR+- These processes will be called 
auxiliary processes. We look for couplings such that the extinc tion probabiliti es for the auxiliary 
processes are easier to compute. For lectures on coupling, see e.g. Lindvalll ( 2002 ). As in Section l3~4"l 
we close the section with a table that sums up our results. 

We denote by (Tk)k^N the successive jump times of the process (TV* , N{ )t>Q, with the convention 
To = and T^+\ = +oo if N^r = N^? = 0. We will be led to consider the continuous time process 
(iV t , iV t 12 )t>Q as well as the discrete time process (N^, N^?)keN- 

Recall also that the time at which the process reaches the horizontal or vertical axes has been denoted 
by r (see 

The denominations for super and subcritical cases are used for small populations but with some 
slight differences with the cases considered in Section even when the reproduction rate is very high 
in comparison of the death rate, there may always be a probabili ty of extinction by dem ographic 
stochasticity. With the terminology of branching processes (e.g. Athreva and Nev . 197dl ). we say 



that we are in the supercritical case if there is a positive probability of survival. In the subcritical 
case, there is almost sure extinction. 



4.1 Wright's model (Model 1) in a small population 
4.1.1 Extinction probabilities: recurrence equation 

A difficulty comes from the fact that the transition rates of (iV t n , A^ 12 )^ g ]R + and the transition proba- 
bilities of the associated discrete time Markov chain vary with the state of the population (see Fig. U]). 
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Figure 4: Evolution of the distylous system (iV/ 1 , -/V 4 12 )i e R + in Wright's model (Model 1). Rates of 
events (for positive i's and j's) are pictured on the left, while transition probabilities of the embedded 
Markov chain are represented on the right. 



Techniques developed in the literature of random walks on positive quadrants usually focus on homo- 
geneo us random walks (e.g. lHaneveld and Pittengerl . 1990; iFavolle et all . Il992l ; iKurkova and Raschell . 
201lh . 



For i,j G N, we denote by pij = Py(3i > 0, N{ = or N^ 2 = 0), where Py means that we start 
with the initial condition Nq 1 = i and Nq 2 = j. By symmetry arguments, we have: 



Vi, j G N, p 



Pj,i 



Moreover, 



when i = or j = 0, p. 



l :J 



1. 



(4.1) 
(4.2) 



We begin with a recurrence equation satisfied by these extinction probabilities: 



Proposition 4.1. (i) The extinction probabilities pij for i,j G N* satisfy the following recurrence 
equations: 



di 



P 



(r + d)(i+ j 
f 



dj 

; ft - W+ (r + d)(i+j) 



Pi,j-1 



2(f + d) 



Pij+l + 



2(f + d) 



Pi+l 



(4.3) 



The family (pij)i,ieN i- s a solution of the Dirichlet problem ft4-3\ ) with boundary condition M-ty . 
Uniqueness of the solution may not hold, but the extinction probabilities (Pi,j)i,jeN define the smallest 
positive solution of this problem. 

(ii) Assume that the probabilities pi t i for i £ N* are given. Then, the other probabilities Pij are 
completely determined. 



P roof. We ref e r to | Lafitte- Godillo n et a j (|20ld ) for the proof, which follows classical arguments found 
Baldi et all (|2002h or iRevuzl (| 19841 ): for (|4.3|) . we consider the embedded Markov chain and apply 



the Markov property at t = 1. ■ 

The recurrence equation (|4.3p may not admit a unique positive solution. Point (ii) of Prop. 14.11 
tells us that to every boundary condition (pi t i, % G N) corresponds a solution. Point (i) tells us that the 
smallest positive solution gives the extinction probabilities pij. The computation of the p ? ; and the 
determ ination of the probabilities p^i's is a work in progress by lLafitte-Godillon. Raschel. and Tran 
(|2010h . 

By coupling and comparison techniques, we manage to define subcritical, critical and supercritical 
regimes, where there is almost sure extinction or positive probability of survival of the random walk 
killed at the boundary. 
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4.1.2 Extinction probabilities: coupling approach 

It is difficult to solve (|4.3p . Our purpose is to approximate the extinction probability with couplings. 
Lower and upper bounds are provided by: 

Proposition 4.2. Assume that (Nq 1 , Nq 2 ) = (i,j) is a pair of positive integers. 

(i) The population goes extinct in finite time almost surely if and only if f < d. 

(ii) Iff>d, there is a strictly positive probability of survival 1 — pij when starting from any i,j £ N* 
and the extinction probability ptj satisfies: 




I 1 1 1 1 1 \ — — i 1 1 1 1 1 \ — 

50 100 150 200 250 300 50 100 150 200 250 300 

Type 1 1 Type 1 1 




50 100 150 200 250 300 0.0 0.5 1.0 1.5 2.0 

Type 1 1 Death rate 

Figure 5: (a)-(c): Simulation of the paths (iV t , iV t 12 )t>o in the critical, subcritical and supercritical 
cases. The number of individuals of genotype {1,1} is in abscissa, the number of individuals of 
genotype {1,2} is in ordinate. For the three simulations, the initial condition is (100,100). (d): 
Simulated extinction probability for (N , N 12 ) in plain line, together with the extinction probability of 
(N U ,N 12 ) in dashed red line and of(N u ,N 12 ) in blue dash-dots. We start with (A^ 11 , iV 12 ) = (1,1). 

To prove Proposition 14.21 we introduce processes that dominate and bound below (iV t n , A^ 12 ) tg K + , 
and that allow to obtain the bounds in ()4.4p . 
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4.1.3 Domination and proof of the lower bound 

In Wright's model, the random walks (|2.24p can be dominated by the random walks N n and iV 12 

f n and Ntf 



with initial conditions Nq 1 and NP, and reproduction rates: 



~ {N} 1 + N} 2 ) instead of T - (iV, 11 + N} 2 ) l N}1>0 . Np>0 (4.5) 

such that until time r, the processes (A^ r, Nj )teR+ an d (A 7 ) } , Nl 2 ) t ^M. + have the same paths (this can 
be obtained by using the same Poisson measures for (A^ 1 , iV t 12 )t e K + and (iV t n , iV t 12 ) te ]R + , see Appendix 
For (Nl\Nl 2 ) tm+ : 

N} 1 =Atf 1 + jf* (L(N? + iVj 2 ) - d x N} 1 ) ds + 

Nl 2 =iV 12 + jf g(A^ + iVj 2 ) - d x iV s 12 ) d S + M] 2 . (4.6) 

We define by N t = N} 1 + AV 2 the total size of the population ruled by As (A^ 11 , A/" t 12 ) t6ffi+ and 

(AT t , N^ 2 )t£R + coincide until time r, they reach the axes at the same stopping time r. Moreover, 
(A^ 11 , A" t 12 )t 6R+ stochastically dominates (A^ 1 , A" t 12 )t eR+ since it allows rebirths of the disappeared 
type on the boundaries {A^ 11 = 0} U {A^ 12 = 0} provided the total population size is strictly positive. 
The coupling presented here also works for general division rates r(N) instead of f. 



The process (N t , N t } )teR+ is a particular case of two-type Galton- Watson process (e.g. lAthreva and Nev 



1970 . Chap. V) where each individual lives during an exponential time of rate r + d. At death, a 



particle is replaced by: 

• zero offspring with probability d/(f + d): this corresponds to the case of a real death. 

• two offspring of the same type as their mother with probability f/(2(f + d)): one is the mother 
and the other is her daughter, of the same type. 

• two offspring of different types {1, 1} and {1, 2} with probability f/(2(f + d)): one is the mother 
and the other is her daughter with the other type. 

Hence, once the trajectories have reached the horizontal axis, for instance, the extinct genotype {1, 2} 
may be regenerated from birth of individuals of genotype {1,2} from individuals of genotype {1, 1}. 
Notice also that the total size of the population N t = N^ 1 + N} 2 is a continuous time birth and death 
process. Let: 

a = inf{t > 0, N t = 0} (4.7) 
be the extinction time of the dominating process. 

Proposition 4.3. Let us consider the processes (AT 11 , N 12 ) starting from the initial condition (Nq,Nq) 
We set N = N^+N^ 2 . 

(i) The total population is a continuous time birth and death process with birth and death rates f and 
d respectively, and: 

E(N t ) = E(J V„)e<-* P(, < +oo | JV = 1) = | (4.8) 

(ii) For the population of type {S 1 , S 2 } E {{1, 1}, {1, 2}} we have: 



E( N rn = — + [E( N r - ^p) ) e -<" (4.$ 
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Proof. Once that it has been noticed that (Nt)t>o is a continuous time birth and death process, 
the results are standar d. It is indeed classical to consider the generating functions for which (e.g. 



Athreva and Nevl . Il97d . Chap. Ill 4-5): 



+00 



F(s,t) = J2HNt = k\N = l)s k 



f(s-l)-e-C F - d )*(fs-<i) r ^ C4 ]_q\ 

s-ft(s-l) - f - _ , V • ^ 

l-fi(s-l) 11 r - «■ 



A;=0 

Since P(cr < +00 | iVo = 1) = lim^+oo F(0, t) and since: 

P(a<i|iV = l)=P(iV t = 0|iV Q = l) = F(0,i)= $ Qzf e ^ )t ) Xf>d 

if f = d 



n 

i+ft 

e (d-f)t-l 

(d-f)t-5 



if r < d, 



we deduce the result by taking the limit in t. The expectations E(A t n ), E(A t 12 ) and E(A t ) are obtained 
by noticing that t 1— > (E(A 4 n ), E(A t 12 )) solves the system (|3.3|) which has been studied in a previous 
part. ■ 

In conclusion, we can distinguish three regimes: 

Proposition 4.4. (i) In the subcritical case f < d, the population (A^ 11 , A r f 12 )t>o goes extinct with 
probability 1 and so does (A/ 1 , Nl 2 ) t >o- 

(ii) In the critical case f = d, the population (Nf , A 7 f 12 )^>o goes extinct with probability 1, and so does 
(Nl , Nf 2 )t>o, but the expectation of the population size E(Aj) remains constant and the extinction 
time is not integrable. 

(Hi) In the supercritical case f > d, there is a positive probability of survival for (A^ , N^ 2 )t^m + , equal 
to (d/f ) N o 1+N o 2 that provides the lower bound in ft4-4\ )- Moreover, in this case, the mean size E(A^) 
tends to infinity with: 

lim ?<M = i. 

t^+00 E(N t ) 2 

Proof. Point (i) is clear. For Point (ii), we use that: 

r+00 r+00 ^ 

E(cr) = / P(cr > t)dt = / = +00. 

Jo Jo 1 + rt 

For Point (iii), we use (|4.8p and the branching property. ■ 



4.1.4 Minoration and proof of the upper bound 

Comparison with an homogeneous random walk To obtain the upper bound in (|4.4p . we have 
to find a process that has an extinction probability higher than (AT t , A t 12 )t G R + . A natural process is 
the following random walk (A t n , Aj 12 )j<=ir + on N 2 : when the process is in (i,j) E (N*) 2 

• individuals of genotype {1, 1} or {1,2} are produced with rate r(i + j)/2, 

• individuals of genotype {1,1} or {1,2} die with rate d(i + j)/2. 

If we consider the associated discrete time Markov chain, we obtain transition rates which are homo- 
geneous: whatever the state £ (N*) 2 , given the occurrence of an event, the chain goes to the 
north or east with probability f/(2(f + d)) and to the south or west with probability d/(2(f + d)) (see 
Fig. [6|). We begin with computing the extinction probability for this process. Then, we show that a 
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Figure 6: Rates of events (left) and transition probabilities (right) of (iV t 11 ' + , A r i 12 ' + ) t£ ]R + and of its 
discrete time skeleton. 



naive coupling between (N{ } , Nj ~ 2 )teR+ an d (N^ , N^ 2 ) te ^ + does not work and exploit the symmetry 
of the problem to obtain the upper bound in (|4.4p . For this, we introduce a coupling of the original 
process (N^ 1 , iV t 12 )t £ R + with two auxiliary processes (N t l,+ , A^ 12 ' + ) 4g K + and (iV* ' , A* 2 ' )teK+- 

Let us give some explanation. (A/ 1 , -A 4 12 ) 4g R + and (A t n , A/ 2 )t e R + have the same birth rates. When 
we are in the upper octant G N 2 , i < j}, then < d(i+j)/2 < dj. In this case, there are higher 

(resp. lower) death rates for the j individuals {1, 1} (resp. the i individuals {1,2}) in (A t , A/ 2 ) te ]R + 
compared with (A* , A^ )teR+- I n the lower octant {i > j}, the reverse holds. Heuristically, for 
the homogeneous random walk (A f , A^/ 2 )t e iR + , the less abundant genotype is pushed more strongly 
towards the axis, which should entail a higher extinction probability. 

Proposition 4.5. For the process (A* , N{ )teR +) we have when r > d: 



dy /dy /d\ i +i 



Hj[ 3tGR + ,Nl 1 = 0orNP = 0) = {^)+{^y-{^) \ (4.11) 



When r < d, this probability equals 1. 



Proof. When r < d, the result is obtain with arguments similar to the ones developed in Section ^. 1.31 
Let us consider the case f > d. We denote by (T^)/^ the jump times of (A^ 11 , A^ 12 ) tGR+ . For the 
associated Markov chain, (iV^)fceN and (N^?)keN move independently. Hence, the survival probability 
is the product of the survival probabilities of each of these processes: 

Py(vt G R+, iV t n > and N} 2 > o) = (l - (l - (J)') . (4.12) 

This hence provides (|4,11[) . □ ■ 

A first naive coupling For the original process (A* , N^ 2 )t^R + as well as for the random walk 
(A* 1 , N^ 2 )t^u + 1 when in the state (i, j), the jump rates to the north and east on the one hand, and 
to the south and west on the other hand, depend only on the sum i + j and are the same for the two 
processes. Thus, a first (naive) idea of coupling of these two processes is to have the jumps to the 
north/east or to the south/west occur simultaneously for both process. All the jumps of the original 
process are copied by the auxiliary process except the following modifications: 

• When we are in the upper octant {i < j}, a jump of (-A t n , A t 12 )i G jR + to the south gives: 

— a jump of (A t , iV t )t e R + to the south with probability (i +j)/(2j), 

— a jump of (A* 1 , Nl 2 )t£R + to the west with probability (j — i)/(2j). 

• When we are in the lower octant {i > j}, a jump of (A t , A 4 12 )i g R + to the west gives: 
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— a jump of (iV/ 1 , A^ 12 )t g R + to the west with probability (i + j)/(2«), 

— a jump of (N} 1 , A^ 12 ) 4e K + to the south with probability (i — j)/(2i). 

One notices that as long as none of the two processes has reached a boundary of the positive quadrant: 
N^+N} 2 = N^+Nl 2 . This entails that the jump rates and times to the north/east and south/west 
remain the same for both processes. 

Unfortunately, this coupling is not sufficient as there exist paths where the original process goes extinct 
before (N t n , N t l2 ) tm+ : see e.g. Fig. [3(a). 



' i 



Figure 7: (a) An example of path for the naive coupling between (Nl 1 , A" t 12 )t g R + (in plain blue 
line) and (iV* , N{ } 2 )teM.+ dashed red line) where the original process reaches the boundary be- 
fore (A 7 "/ 1 , iV t . The two processes are started in the lower octant {i > j}. For the first step, 
the west move of the original process is replaced by a south move for the homogeneous random walk. 
Then, the original process keeps moving west, and since it is in the upper octant {i > j}, these moves 
are the same for the auxiliary process. The original process reaches the vertical axis first, (b) When 
we use the coupling with the two processes {N t 1,+ , N t 2 ' + )tem + and (A r 4 11 '~, A r 4 12 '~) t6 K + (dashed red line 
started in the upper and lower octants respectively), things happen as follows. For the first move, 
since the original process is in the lower octant, it is coupled with (N t l '~ , N t 2 '~)t£M. + and the west 
move of the original process is here replaced with a south move of the latter process ; by symmetry, 
the auxiliary process (JVj 1 ,+ , A^ 12 ' + )£ 6 ]g + in the upper octant makes a jump to the west. After this, 
the process (iV/ 1 , AT t 12 ) te jR evolves in the upper octant and is coupled with (iV t 11,+ , N^ 2 ' + )f^g,: all the 
west moves are unchanged. The auxiliary processes reach the boundaries first. 



Coupling with two copies of (A 7 ^ p , N{ } 2 )teM.+ To prove the upper bound in (|4.4p . we introduce a 
coupling of (Nl\Nl 2 ) tm+ with2 copies (7V t n ' + , A> 12 '+) t6R+ and (ivf - A> 12 '") ie iR + of (A> n , A> 12 ) ieM+ . 
Their respective initial conditions are: 

(iV U ' + ,iV 12 ' + ) = (minK 1 ,^),^^ 11 ,^ 12 )) 
(iV n '-,iV 12 '-) = (m^lJVo 11 ,^),!!!^ 1 ,^ 2 )). 

As for the process (A"/ 1 , N^ 2 )t^u.+ , there exists a coupling (see Appendix I A. 2 j) such that: 

• the two processes (A^ 11 '" 1 ", A^ 1 ' + )t£M. + and (N t , N t 2 '~)teR + are symmetric with respect to the 
line {i = j}, and started respectively in the upper and lower octants, {i < j} and {i > j}. 
Thus, there is at every time one of these processes in each octant {i < j} and {i > j}. Of 
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course, since these processes are copies of (iV/ 1 , iV t 12 )t e K + they may change octant, but this is 
done simultaneously and they meet on the line {i = j}. Moreover, these processes both reach 
the axes at the same time. 

• the jump times to the north/east or south/west are the same for (JVj , A r i 12 )j g K + and the two 
auxiliary processes. Indeed, when in the state these jump rates are f(i + j) and d(i + j), 
and as long as the north/east and south/west jumps of the three processes occur simultaneously: 

A^ 11 + N} 2 = N^' + + Nl %+ = Nt 1 ' + Nl 2 '~. 

• the coupling is between the original process (iV* , N{ )teM.+ an d the auxiliary process that belongs 
to the same octant, and it depends on this octant. The other auxiliary process is obtained by 
symmetry. If at a time t, A^ 11 < N^ 2 (we are in the octant {i < j}), the coupling determines 
as follows the behavior of the auxiliary process that is found in this octant {i < j} at t (when 
N} 1 = Nl 2 , we choose (N t 11,+ , jV t 12,+ ) by convention): 

— If the original process moves to the north, east or west, so does the chosen auxiliary process. 
If the original process moves to the south, then the chosen auxiliary process moves to 
the south with probability (A^ 11 + N} 2 )/{2Nl 2 ) and to the west with probability (iV t 12 - 
Nl l )/{2Nl 2 ). 

— The behavior of the other auxiliary process is obtained by symmetry with the first bisector. 

If at time t, iV"/ 1 > A^ 12 , the coupling determines the behavior of the auxiliary process that is 
found in this octant {i > j} at t as follows: 

— If the original process moves to the north, east or south, so does the chosen auxiliary process. 
If the original process moves to the west, then the chosen auxiliary process moves to the west 
with probability (iV t 11 +^ 12 )/(2^ 11 ) and to the south with probability (Np - N^ 2 ) / (2Np) . 

— The behavior of the other auxiliary process is obtained by symmetry with the first bisector. 

It can be easily checked that the auxiliary processes have the same distribution as (A^ 11 , A^ 12 )t G R + 
(see Appendix) and that for every time t, as long as neither of the three process has reached the 
boundaries {i = 0} or {j = 0}: 

min(A> f n ' + , Nl 1 '-) < N} 1 < max(A\ 11,+ , jvj 1 ") 
min(A> f 12,+ , Ml 2 '') < Nf 2 < max(A> t 12 ' + , N^''). 

As a consequence, 

{3t G R+, Np = or A^ 12 = 0} 

C{3t £ K+, iV f 11,+ = or N~t %+ = 0} U {3t G M + , A^ n '~ = or n] 2 ^ = 0} 

={3t G M+, N^' + = or = 0}, (4.13) 

since the two processes (A^ 11 '" 1 ", A" f 12 ' + )t g iR + and (A^ 11 ' - , A" t 12 '~)t g iR + reach the axes at the same time. 
The computation of the probability of the event on the r.h.s. of (|4.13p is given by Prop. 14.51 and 
proves the upper bound in (|4.4p . 

4.2 Fecundity selection model (Model 3) in a small population 

For the fecundity selection model, we will prove that the same criteria as in large population hold for 
separating the subcritical, critical and supercritical cases: 
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Figure 8: Evolution of the distylous system {N} , N t )teR+ i n the fecundity selection model. 



Proposition 4.6. Assume that (Nq 1 ,^ 2 ) = (i,j). 

(i) If r < 2d, then we have almost sure extinction in finite time. 

(ii) If r > 2d, then we have a positive probability of survival 1 — pij such that: 



(- 



d\i+j /2d\ min (*.i) 

<Pi,j< — ■ (4-14) 



Proof. We begin with the proof of (i). We first assume that r < 2d. When in the state £ 
N 2 \ {(0, 0)}, with i < j, the total death rate of the population is d(i + j) and the total birth rate is: 

ij j i _/l \j — i\ \ _/l \j — i\ \ 
2r- :=r- : i + r- :J =r{ — I — ^") z + r ( 



i + j i + j i + j J v 2 2(i + jY v 2 2{i + j) JJ 

f , . r 1 7 j) f , 

=-(* + J + ~ — < -(i + j)- 

2 V JJ 2 i + j - 2 V Jy 

Hence, it is possible to dominate the total population size by a continuous time Galton- Watson process 
with individual birth and death rates f/2 and d. For f < 2d, the latter process is critical or subcritical 
and extinction is almost sure. 

We now prove (ii), and assume that f > 2d. Since the birth rates are bounded above by fi and 
fj, (iV/ 1 , Nt )t£R + is bounded above by the process (N^ , N t ) of Section T4.1.3I We turn to the 
minoration. When (iV t , A^ t 12 ) = (i,j) with i < j, then the birth rate for the population {1, 1} is: 

- ij . f . 



i + j- 2 

Hence, as long as N\ < A^ 12 , A^ 11 is stochastically lower bounded by a supercritical continuous birth 
and death process with individual birth and death rates f and d. Similarly when j < i, N 12 is 
stochastically lower bounded by the same process. 

As a consequence, the total size of the population (A^ 11 + A" t 12 )i e K + is stochastically lower bounded by 
the supercritical continuous time birth and death process introduced above, started at min(A r Q 1 , A^ 2 ), 
and which survives with probability 1 — (2d/f) min ( N ° ' N ° ' (see Prop. 



Remark 4.7. Notice that the upper and lower bounds in \J^.lJ\ are not tight and the couplings used in the 
proof are not likely to be optimal. For the lower bound, the bound 1 for i/(i + j) and j / (i + j) is rough. 
For the upper bound, we have only considered the less represented component in the population. ■ 

4.3 Compatible population without pollen limitation nor inbreeding depression 
(Model 4) 

As a reference for the previous simulations, we provide the results for populations where there is no 
self-incompatibility. This corresponds to the case where Vx 6 W 1 , <£(x) = 0. In this situation, we 
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can forget the phenotype of the individuals. This leads us to consider a population where individuals 
reproduce with an individual rate R(N) and die with a rate d, N being the size of the population. 
Notice that this corresponds to the process (N t )t>o studied in Proposition 14.31 

The Markov chain embedded in the continuous time branching process is the chain on N, with 
transition probabilities qij from i to j defined for i > by: 

= nfl'l 3 > = p /.f , 90,0 = 1 and else = 0. (4.15) 

it(«j + d R(i) + d 

Let ^ be the extinction probability when the initial population is of size i. 

Proposition 4.8. Let pi be given. Fori > 1, 



Pi+i = Pi 




Proof. Using the strong Markov property at the time of the first event (see e.g. iBaldi et all . 120021 ) : 

d R(i) ,~ 

Pi = RMTd Pi - i+ mTd Pi ^ (417) 

where by convention po = 1. We deduce from this that: 

Pi+l ~ Pi =J^(Pi-Pi-x) = n / fl(j) frl " !)" 

by recursion. The result follows by using = p\ + 5Zj=i(Pj+l ~ Pi)- ' 
Example 4.9. In the case where the individual reproduction rate is constant R(i) = f, 



In this case, (Nt) t >o is moreover a one- dimensional continuous time Markov branching process (e.g. 



Athreva and Neu . \197u . Chap. Ill) and we know that p\ solves g(s) = s where g(s) is the generating 



function of the offspring distribution: 



9(s) = -r— + 



d + r d + r 
This gives 

d ( d\ ' 

Pl = — and hence Vi > 1, Pi = f — J , (4.20) 

which is expected since the branching property holds in this case, contrary to cases where R(.) is not 
constant and where there is interaction between the individuals. 

In a nutshell, we have obtained the following estimates for the Pij's: 
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Case 


Case description 


Extinction probabilities (or bounds) when started from (i, j) 


Wright's model (Model 1) 


(a-b.l) 
(c.l) 


f <d 
r > d 


Pi,j = 1 

{djr) i+ i < PiJ < (d/f) 1 + (d/ry - (d/f) i+ i 


Fecundity selection model (Model 3) 


(a-b.3) 
(c.3) 


f < 2d 
f>2d 


Pi,j = 1 

(d/f)*+J < Pij < (2d/f) min( - i ^ 


Self-compatib. 


e model without inbreeding depression (Model 4) 


(a-b.4) 
(c.4) 


r < d 
r > d 


Pi,j = 1 

PiJ = (d/rf+K 



Table 2: Summary of the bounds on the extinction probabilities pij defined in Section \4-l.l\ It is seen that 
the disjunction between sub, super and critical regimes is the same as for large population. In the super- critical 
case, there remains however always a positive probability of extinction. 



5 Simulations 

5.1 Simulation algorithm 

The population dynami cs described in the previous Sections 12.11 and 12.21 can be simulated with the 
following algorithm (see Fournier and Meleard . 20041 ). Notice that the algorithm that we propose is 



exact (in the sense that it describes exactly the dynamics described above without approximation 
scheme). Assume that the population is known at time t. Then: 

1. We define the total event rate at the population level by: 

C t = R(N? V ,Nt)N? + N t d. (5.1) 

l<u<v<n 



2. The next event time is if = t + r where r is an independent random variable that is exponentially 
distributed with parameter Ct- 

3. We then draw an independent uniform random variable 9. 

If < 9 < Ei< u <v<n R(N? V ,N t )Nt V /C t then a birth happens: 

(a) The ovule is of type {u, v} with probability 

R(N^ v ,N t )N^ v / Yl R(N^ V ,N t )N^ v . 

l<u<v<n 

(b) The pollen is then of type {u',v'} with probability p% v (u',v'). 

(c) The offspring is then of genotype {n, u'}, {u,v'}, {v,u'} or {v,v'} with probability 1/4. 

If 9 > J2i< u <v<n R\Nt i Nt)N t /Ct then an individual dies. This individual is drawn uniformly 
among the living individuals. 



5.2 Simulations performed in the case of a distylous species 

Each simulated curve is obtained as the average on 5000 simulations of paths (iV t , A r i 12 ) tg [ 0)10 ooo] • The 
estimated extinction probability is obtained as the mean number of extinctions before time 10000. The 
simulations were run for varying number of different genotypes as initial conditions when considering 
Wright's model (Model 1, Section [2. 2p . while Nq 1 = Nq 2 = 1 when considering the dependence model 
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(Model 2). 



In Fig. [9] the extinction probability for a constant reproduction rate r = 2 is shown for an in- 
creasing death rate, for both the distylous and the self-compatible case, for different initial numbers 
of the two genotypes, under the Wright's model. Fig. [9] shows that in small population, the extinction 
probabilities are higher for the distylous than for the self-compatible population. As expected, the 
difference between the distylous and the self-compatible populations is lower when the initial popu- 
lation size increases, as shown in proposition 13.11 for Wright's model. The results shown for Wright's 
model shows the importance of the absorbing effect that increases the extinction probability for the 
distylous population. We also see that when a < +oo in (|2.7p . the extinction probabilities are higher 
than in Wright's model, both for distylous and self-compatible populations. This reflects that when 
pollen limitation is high, it is more difficult to encounter a mate and to produce offspring. 

In Fig. [TOl we estimated the ratio p of the extinction probability in the self-compatible case on the 
extinction probability in the SI case, which are obtained from our simulations results. The Fig. [Tola) 
shows that the ratio p gets smaller when the initial population size is higher, which means that the 
extinction probabilities for the self-compatible population decreases more rapidly when N increases 
than for the distylous population. This effect is even worse when the initial population is asymmetric, 
that is when a given genotype is more frequent than the other (compare filled with empty symbols). 
Once again, the differences between distylous and self-compatible populations shown in Fig. [TUTa) are 
only due to the absorbing effects. Fig. fTOTb) shows the ratio p for different initial conditions under 
the dependence model with a = 1 and j3 = 100. Fig. [TUTcl shows similar results but for different 
strength of pollen limitation (a < +oo) with iVg 1 = ./Vq 2 = 1. It is remarkable in this figure that 
when the strength of the pollen limitation is high (a is small and /3 is large) then the ratio p is higher, 
which means that the higher the pollen limitation, the lower the difference between distylous and 
self-compatible populations. 

6 Discussion 

Impact of pollen limitation versus demographic stochasticity and boundary effects Thanks 
to the three relationships between the compatible population size and the reproductive rates we as- 
sumed (see 12.21 and Fig. [1]) , we are able to disentangle the relative effects of demographic stochasticity 
and pollen limitation on the fate of the populations. Indeed, when Wright's model is assumed, there 
is no pollen limitation. In other words all individuals of a given mating types receive enough pollen 
to fertilize all their ovules, as long as at least one compatible individual is present in the population. 
This assumption introduces a discontinuity in the individual rate of seeds production. On the other 
extreme, the fecundity selection model has a rate of seeds production that vanishes continuously at 
the boundaries. Finally, we investigated the combined effect of pollen limitation and demographic 
stochasticity thanks to the dependence model. 

We have first studied large populations. We exhibited three different regimes, depending on the 
relationships between the birth and death parameters. Under the subcritical regime, the population 
goes extinct when t — > +oo, while its size tends to infinity under the supercritical regime. As expected, 
when there is no pollen limitation, a distylous population behaves like a self-compatible population. 
When there exists an equilibrium, the population size is n* in both cases. To the contrary, under the 
fecundity selection model, the size at equilibrium is In* showing that pollen limitation can have a 
large effect on the dynamics of distylous populations, even in large ones. In the dependence model, a 
saddle point appears that makes condition for the maintenance of the population more scarce than in 
Wright's model and favors populations with symmetry in the initial conditions. 

In small populations, the separation into subcritical, critical or supercritical regimes is the same 
as for large populations. In supercritical cases, we are able to find lower and upper bounds of the 
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(a) (b) 
N t1 = 1,N 12 =1 Nn=2,N 12 =2 
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Death rate, d Death rate, d 

(c) (d) 
N 11 = 5,N, 2 = 5 Nn = 10, N 12 = 10 




0.5 1.0 1.5 2.0 0.5 1.0 1.5 2.0 



Death rate, d Death rate, d 

Figure 9: Comparison of the extinction probabilities of (Nj Nl 2 ) t >o for a constant reproduction rate 
f = 2 when the initial number of the different genotypes varies (indicated in the top of the box). Thick 
line: Extinction probability in the self- compatible case (from equation \4-20\ ; Crosses: self- compatible 
case with pollen limitation (a = j3 = 1) ; Full circles: distylous case with Wright's model ; Circles: 
distylous case with the dependence model (a = (3 = 1). 



probability of extinction of a distylous population under the Wright's and fecundity selection models. 
Again, the main difference between these models relies on the fact the critical condition for the survival 
of the population is twofold higher in the fecundity selection model. 

Our simulations show that when pollen limitation is low (a = 1 and (3 = 1, see Fig. [9|), the increase 
in the probability of extinction is small relative to the impact of demographic stochasticity. When 
we compare the extinction probabilities between a distylous and a self-compatible population by the 
measure p, we see that the higher the pollen limitation (when a decreases), the lower the difference 
between them since p increases (see Fig. fTOj) . Those results have important ecological and evolutionary 
consequences since they suggest pollen limitation plays a minor role in small populations relative to de- 
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Figure 10: Ratio of extinction probabilities in self- compatible cases over distylous case with r = 2, 
varying initial conditions and varying a and f3. (a): Wright's model, Nq 1 = Nq 2 = 1 (triangles); 
Nq 1 = iV 12 = 2 (empty squares); iV n = iV 12 = 5 (empty circles); iV n = 3 and iV 12 = 1 (filled 
squares); Nq 1 = 9 and Nq 2 = 1 (filled circles), (b): Dependence model with a = 1 and (3 = 100 and 
varying initial conditions: N^ 1 = Nq 2 = 1 (triangles), Nq 1 = Nq 2 = 2 (squares), Nq 1 = Nq 2 = 5 
(circles), Nq 1 = Nq 2 = 10 (crosses), (c): Dependence model with Nq 1 = Nq 2 = 1 and (3 = 1 (filled 
symbols), f3 = 100 (empty symbols), a = 0.1 (circles), a = 1 (triangles), a = 2 (squares). 



mographic stochasticity, and especially the stochastic loss of one of the two mating types. When pollen 
limitation is very large (when a is low), small populations of self-compatible and self- incompatible 
species tend to behave similarly. In short, our results suggest that pollen limitation plays a major role 
in large population only. 

Mate Finding Allee effect, inbreeding depression and the evolution of SI One of the 

most intriguing and long-standing problem in evolutionary biology resides in the existence and main- 
tenance of SI. Indeed, the establishment probability is lower in the case of SI compared to self- 
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compatible species s ince a single individual is sufficient to colonize empty spaces (the Baker's law, 
Pannell and Barettl (|1998|)). Furthermore, as shown here, the extinction probability is higher for SI 
than for self-compatible species in absence of inbreeding depression. It has been also shown that in 
infinite populations, the conditions for the invasion of SI po pulations by a self-compati ble mutant 
are less stringent when there is strong pollen limitation (e.g. Porcher and Lande . 2005bl ). The pro- 
cess generally invoked to explain the existence and maintenance of SI is inbreeding depression: SI 
can be advantageous relatively to self-compatibility when the cost of inbreeding depression caused 
by self-fertilization is high. Inbree ding depression can play a role in two ways: by preventing the 
invasion of self-compatible mutant ( Porcher and Lande. 12005b), or by increasing the extinction rate 



1995), what is suggested in data since 



of obligate selfing species relatively to SI species (jLvnch et al 
self-compatible species are most often localised at the leaves of ph ylogenies, esp ecially in plant fami- 
lies w here distyly is prese nt (e.g. in Amsinckia ISchoen et all (|1997l ). in Narcissus iPerez-Barrales et al 



and in Psvchotria ISakai and Wrightl (|2008l )). Here, we proposed a very simple model taking 
into account inbreeding depression caused by self-fertilization, to investigate if there are conditions 
under which the extinction probabilities are higher for SC populations than for SI populations. We 
found that there are some conditions in large populations where the size at equilibrium is lower in SI 
populations than in self-compatible populations: under strong inbreeding depression, SI populations 
may be less sensitive to extinction. 



What about more complex dominance networks? Although we developed a general model for 
the population dynamics of SSI species, we mainly investigated the dynamics of a distylous species. 
Distyly is however the case where the impact of the existence of mating types is the highest since 
only two mating types exist. When the number of mating types increases, the proportion of com- 
patible individuals also increases. It would be interesting to investigate the impact of the dominance 
relationships among S-alleles on the extinction of population, to check if there are dominance in- 
teract ions patterns that are less sensitive to extinction than others, as highlighted by Kirchner et all 
(2006): :hey showed that populations where all S-alleles are codominant have a higher extinction rate 
than po pulations wh e re a l inear hierarchy of dominance exist between S alleles (the DOM model, 
see e.g. (Billiard et a3 d2007f >). Our model could be used to investigate the effect of these dominance 
interactions more precisely. 



A SDEs and Proofs of the Section 12.31 Propositions 
A.l SDEs 

Following IFournier and MeleardI <|2004h . we present a SDE describing the evolution of (Zt)t€U. + - 

Definition A.l. Let Q(ds,dg' ,d9) be a Poisson point measure on x E x M + with intensity 
q(ds,dg' ,d0) = dsdn(g')d9, where ds and d9 are Lebesgue measures on E + and where dn(g') is the 
counting measure on E. To each atom from Q(ds,dg' ,d0) are hence associated a time of possible event 
s, the genotype g' that either appears or dies and an auxiliary variable 9 that decides what happens 
(with a role similar to the variable 9 in Point 3 of the algorithm of Section \5. 



Z t {dg) = Z (dg) + I [ 6 g ,(dg) U e<rg > (Zs 

JO JExM.+ K ~ 



and where r 9 ' (Zt) has been defined in h2.10]) . 



.) 
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rs' (Z s _ )<6<rs' (Z s _ )+dxNj>_ 



Q(ds,dg',d9), (A.l) 



Existence and uniqueness of a solution of (|A.lj) are stated in the next proposition. Moreover, it is 
possible, for a given test function /, to derive from (|A.ip equations for the evolution of ({Zt, f))t 
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Proposition A. 2. (i) 7/E(Ao) < +00, then there exists a unique solution to SDE IIA.1)) . 

(ii) If additionally E(Aq) < +00 then for any bounded test function f on E, 

(Z t ,f)=(Z ,f)+ f Y, {r uv (Z s )-dxNDf{{uM)ds + M[ (A.2) 

JO r..3t^ 



{u,v}£E 

rf 



where (M/)j<=r is a square integrable martingale starting from with quadratic variation: 

(Mf) t =f t Y {r u "(Z s )+dxNr)f 2 (u,v)ds. (A.3) 
{u,v}eE 

Proof. The rate r uv defined in (|2,l(jp is bounded by a linear function in Aj. The moment condition 
in (i) allows us to prove tha t there is no explosion. The proofs then follows the ones developed in 



Fournier and Meleardl (|2004l . Th. 3.1 and Prop. 3.4) for a model of plant with asexual reproduction. 



A.2 Couplings in Wright's model for the proofs of Section 3] 

In this section, we give the expression of the processes ( A t n , N^ 2 ) t £R + , ( A t n , A t 12 ) tg n + and ( A t n , A t 12 )t G R + 
that appear in Section HI 

Let us rewrite Equations (|2.24p thanks to (jA.lj) : 

= N U + / / V={1,1} ( >0;iV s 12 >ol f(iv]i +iv]2 ) 

Jo JexR+ v " e< . 



- IfwJi +nP ) f(jvJi +ivp ) ) Q(ds, dg', dO) 

=2 —<S< =5 —+dxN}l ' 

and similarly for (A/ 2 )j G ]r + . 

For the stochastic domination, the process (A/ 1 , A/ 2 )t e K + introduced in (|4.5p can be rewritten as: 

Nt 1 = No 1 + f [^ m \'={i,i}{^ n ^f(N}i+Np_) 

)Q(ds,dg',dO) 



^5 — <9< == —+dxNl 1 



and similarly for (A t 12 )t g R + . If we start at a point where (A/ 1 , A f 12 ) = (A/ 1 , A/ 2 ) then it is clear that 
these processes have the same births and deaths. If one of the components is null, then there is no 
birth in (A/ 1 , A t 12 ) ieM+ while there is still births in (A t n , iV t 12 ) t6R+ . 

For the stochastic minoration, we have used in Section T4. 1.41 the processes (A^ 1 '" 1 ", A 4 12 ' + )t e ig + and 
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(A f n ' , A f 12 ' )teK + . Let us give an SDE to describe their dynamics. 



\n}1 -N}1 )(7V"' + -N™> + )>0 X 



A^ + =min(iV ^iV 12 ) + / / 

JO JEx] 

{ljVii<JVi2 ^%ii_jv s 11 >o - % 8 U -JVJ-1<0 

- V={l,2}l(r+d)(iVll +NP ) rtJVs 11 +iVl 2 ) ) 

f — <#< =2 — / 

+ fl/V" >Af s 12 ( IjVj-l-iV" >0 - ^g'={l,l}^f(NP_+NP ) (r+d)(JVli+JVl£ ) ) \ 
— — \ — ^fl< / J 



"a <^ 5" 



AT11=AT12 ^ATii_7Vii >0 - J1AT11-AT11 < 
"^(Am -JVJ-2 )(7V s 1 ^' + -7V s 1 ^+)<0 X 

{^11 <iVi2 (ljV-i2_jvi2 >0 - V={1.2}^« 



-<6»<- 



+ ljV,u >Afi2 I ljvja_m2 >0 - %i2_jvja 



<0 



2" — <#< -2 —+dN}i ' > 

Q(ds,dg',d9) (A.4) 

Similar SDEs can be written for (N t 2 ' + )teR + , (^Vf 11 ' _ )tgiR + and (Nl 2 '~) t ^m + - Notice that the process 
(A"/ 1 , NP) t £M + appears in the definition of the auxiliary processes because of the coupling. 
The behavior of (A t ' )teR + depend on whether (A f n , A f 12 ) fe R + and (A t 11,+ , A^ 12 ' + ) tg K + belong to the 
same octant and second on which octant it is. These processes belong to the same octant at time t if 
(A"/ 1 - A t 12 )(A t 11,+ - Nl 2 ' + ) > 0. In this case (3 first lines of {53$): 

• Births of individuals {1, 1} in the auxiliary process occur as soon as an individual {1, 1} is born 
in the original process (iV t , A t 12 )t g K + , 

• If we are in the octant {i < j}, then every death of individuals {1, 1} for the original process 
entails a death of an individual {1,1} for (N t ' + , A* 2,+ )t 6 R + . Additionally, some deaths of 
individuals {1,2} for the original process are changed into deaths of individuals {1,1}. This 
happens (second line of (|A.4j) ) when: 

(f + d)(N£ + N™_) n ^ r(J\£ + N?_) 



2 

which corresponds to a rate 



< 9 < ~ — + dN 



f(N^+Nj 2 _) ^ _ (r + d)(N^+N^_) = d(N}*-N}*) 
2 s - 2 2 

Hence individuals {1, 1} in iV 11,+ die with rate: 

d(N?_ - N^_) d(N?_ + A s 12 ) d(N}b + + Al 2 ' + ) 

2 s ~ 2 2 

In this case, the process (N t ' + , N t 2 ' + )teR + has the rates announced in Section [4.1.41 

• If we are in the octant {i > j}, then deaths of individuals {1, 1} for (N t ,+ , A t 12 ' + )f g K + originates 
from the deaths of the original process such that: 



r(N^_ + At 2 ) (f + d) (Nil + Nj 2 ) 

< 9 < . 

2 ~ 2 
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The death rate for the population {1, 1} is thus d(N}_ + N]_)/2. 

When the original process (A^ 11 , A^ 12 )f g K + is on the line {i = j'}, the death rates are d = d N^ 2 = 
d(N^ +N} 2 )/2 and no correction is needed for (JV t n ' + , A^ 12 ' + )i eR+ . 

Finally, when the original process (A^ 11 , AT t 12 ) t6R+ and the auxiliary process (N t ' + , A^ 12,+ )t 6 K + do not 
belong to the same octant (lines 7-9 of (|A.4[> ). the auxiliary process (N t '~ , N t 2 '~)t£M.+ is constructed 
as above and the process (A^ 11,+ , A^ 12,+ ) tgR+ is deduced by symmetry. 

A. 3 Large population limits and proof of Proposition 12.21 

Sketch of Proof of Prop. \2.Si We begin with (ii) by assuming existence of a solution to: 

&,/> = <&,/>+ f\ E r uv (Z s )f({u,v})-d(Z s ,f)\ds, (A.5) 

\{«,»}6E / 

where / is a bounded test function on E. Since all finite measures of Mp(E) have the form (|2.19p . it 
remains to prove that the functions n uv satisfy (|2.20p . Let us choose f({u,v}) = lr u i tV n({u, v}) and 
integrate / with respect to Since (£t, Iwyi) = nf' v> , (|A.5[> gives: 



+ / {r u ' v '(Cs)-dn^ v ')ds. (A.6) 



Notice that for any finite measure ( on £ of the form (|2.19p . r uv (£) is the sum of terms of the form 
fn uu p uu (v, v') for u, u' , v, v' in [1, nj. In Wright's model, there may be discontinuities when several of 
the possible genotypes reach a size 0. Else, the function r uv (.) that we consider is locally Lipschitz con- 
tinuous as products of locally Lipschitz continuous functions. Since £ is continuous and since r uv (.) 
is continuous, the integrand in ()A.6p is continuous, which entails that n u v is of class C 1 and hence 
of class C°° by direct recursion. Taking the derivative with respect to time gives (|2.20p and achieves 
the proof. Moreover, the conditions on r uv imply by the Cauchy-Lipschitz theorem that there exists a 
unique solution to (|2.20p and hence to (1A.5[) . Let us now prove existence, which is a consequence of (i). 

First, let us notice that for given K G N* and real test function / on E, the process (Z^)teR+ 
satisfies the following evolution equation: 

(zW,f)={4* ) ,f)+ f[ £ r™{Zf))f{{uM)-d{Zf\f)\ds + M^ (A.7) 

J ° \{u,v}GE j 

where (M^'^t^M. is a square integrable martingale starting from with quadratic variation: 

<MW> t =if E r™(zW)f({u,v}) + d(Zf\f) ) ds. (A.8) 

Heuristically, as the quadratic variation of the martingale is of order 1/K, the stochastic part of the 
process will disappear in the limit. Since the jumps of are of order 1/K, the limiting values of 
(ZW)xeN* are necessarily continuous. Moreover, (|2,18p implies that every one-dimensional marginal 
has a finite mass. The limiting values hence belong to C(M + , Aip(E)). 

The proof of (i) separates classically in two steps. We establish tightness of the l aws of (Z^)kpn 



which implies that this family of probability measures is relatively compact (e.g. lEthier and Kurtz 



19861 . p. 104). Then, we establish that every limiting value solves (|A.5P which has a unique solution. 
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For any T > 0, the tightness on D([0, T], M.f{E)) is obtained by using a criterion due to iRoellv 



( 19861 ). and since E is finite, the problem amounts to prove the tightness on B([0, T], M) of the sequence 



((Z( K \ f))K£N* for bounded test functions /. Given the local Lipschitz continuity of r uv (.), given that 
|r u1, (Z)| < CfNt and given the mo ment estimates (I2.18p thi s is a clas sical computation which uses 



Aldous and Rebolledo criteria (e.g. Joffe and Metivierl . 19861 ) . See e.g. Fournier and Meleard (2004, 
Proof of Th. 5.3). 

Taking the limit in (|A.7p thanks to (|2.18p again allows us to identify the adherence values of 
(Z( K ')x;eN* as the solution of (|A.5p . This provides existence of a solution to (|A.5|) and since we have 
established uniqueness, there is a unique limiting value to which the sequence converges. ■ 
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